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On Modified Divided Differences I 


Although divided differences of a function are of basic importance in the 
theory of numerical analysis, they are not nearly as useful as ordinary 
differences in application. In the first place, they are much harder to gen- 
erate; not only because they involve unequal intervals of the argument, but 
because divisions are involved, often by small numbers. Furthermore, they 
are harder to interpret than ordinary differences, and the problem of making 
judgments at each step about the number of meaningful significant figures 
in the successive differences is one of considerable difficulty for high-speed 
computing machines, or for I.B.M. machines. We shall show that part of 
this difficulty at least can be overcome by modifying the definition of 
divided differences. These modified divided differences take on many of the 
characteristics of ordinary differences, based on equal intervals between 
successive arguments, and the familiar theory relating to ordinary differ- 
ences becomes a special case of the more general theory. Moreover, the 
decimal point can remain fixed, and all differences can be carried to the 
same number of decimals as the entries themselves. Such modified differ- 
ences can therefore serve as useful tools in analyzing functions which are 
available only at unequal intervals of the argument. 

Definitions. Let a function u(x) be given for a sequence of arguments 
Xo, X1, ***, Xn, all included in a closed region X. No assumptions need be 
made about the ordering in the sequence. For brevity, let u = u(x,). The 
accepted definition! for the first divided difference of u, and 1; is 


Un — Uj 

[xx;] = 

The divided difference of order m, involving uo, %1, U2, ---, Um is obtained 
from the divided differences of order m — 1 by the formula 


It is well known that [xox; --- Xm] is a symmetric function of the arguments. 
If u(x) possesses m continuous derivatives in X, then 


(1.1) , tink. 
In particular, if u(x) is a polynomial of degree m, then the mth divided 
difference is a constant, given by the right-hand side of (1.1). In the follow- 
ing discussion it will always be assumed that the function under considera- 
tion has ” continuous derivatives in X, whenever the nth difference (ordinary 
or divided) is referred to. 

It is to be noted from (1.1) that the order of magnitude? of a divided 
difference does not depend strongly on the size of the intervals between 
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2 ON MODIFIED DIVIDED DIFFERENCES I 


successive arguments. In contrast to this, the mth ordinary difference Au 
is related to the mth derivative of the function as follows: 


(1.2) A"uy = sin kX, 

where hk is the constant difference between successive arguments. Hence 
| A"uo| decreases with h; this is a desirable property. In practice, # can often 
be chosen small enough so that an approximating polynomial of a given 
degree shall represent u(x) with desired accuracy. We can introduce a 
modified divided difference possessing similar characteristics as follows: 


Define the modified first divided difference of uo and u; by 
w 
(1.3) [w; ] m wl = (uy — Uo), 


where w is an arbitrarily chosen positive constant. We shall show in Section 
II that if there are NV arguments x; in X, then a useful constant to choose 
is one which is close in magnitude to 


If the arguments x, form an increasing or decreasing sequence in X, then 
(1.4) becomes 


(1.5) w= |xv-1 — x0|/(N — 1). 


In many cases the order of magnitude of x, — x,-1 remains the same over 
the region X, even though the intervals are not all equal. If w is chosen as 
the average length of the grid intervals in X, w/(x; — xo) will usually be 
of the order of magnitude of unity. In that case the number of decimals in 
the first modified divided difference which have meaning is the same as in 
the entries u(x) themselves. 

We define 


{[w; +++ — [w; Xn 


(1.6) [w; xox. xn] = 


Again nw/(x, — xo) will usually be of the order of magnitude of unity, and 
it is reasonable to retain the same number of decimals in all the divided 
differences as there are in u(x). 

Since w was introduced as a device to permit the retention of a fixed 
decimal point, it is enough to use a value of w which is close to, but not 
necessarily equal to, that defined in (1.4). 

It follows from (1.1) and (1.6) that 


(1.7) [w; +++ Xn] = tin X. 


for 


(1.8 
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ON MODIFIED DIVIDED DIFFERENCES I 3 


The character of the modified divided differences parallels that of ordinary 
differences, and the rate at which the successive modified divided differ- 
ences fall off again permits an estimation of the accuracy of a given approxi- 
mation formula for u(x) over X. Let 


x — Xo = pw; Xe — Xo = 


The well known identity for divided differences now assumes the following 
form : 


= wo + plw; + PY 
+ p(p — [w; + 
(1.8) 
pi)(p — -(p a8 [w; xn] 


If the right-hand side of (1.8) is truncated after the divided difference of 
order m, then the remainder R can be expressed by 


d**1y (t) 


(1.9) R = M[w; xxox, --- x, ] = Mw aun’ 


where 


tin X, 


M= +1)! 


Whenever M is numerically less than unity, the magnitude of the divided 
differences of order (m + 1) provides an estimate for an upper bound of the 
truncation error R. 

We can write 


(1.10) [w; Xml = Ds 
Il (xz — x5) 
where the prime indicates that the factor corresponding to j = k is to be 


omitted. In the special case of a divided difference of even order 2n, chosen 
centrally around xo, (1.10) becomes 


wm 


(1.12) 
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In the case where all intervals are equal to w, p, = k, and (1.12) reduces to 
the ordinary central difference. 

Since the only difference between the modified and ordinary divided 
difference of order m is the factor wm! for all entries in a column of differ- 
ences, it would appear that just as much information could be obtained 
from the entries without modification, provided the proper number of sig- 
nificant figures were retained in each difference. This is true. The chief 
advantages of modified divided differences over ordinary divided differences 
can be summarized as follows: 


1) If wis suitably chosen the modified differences can be retained to the 
same number of places as the entries u, themselves. 

2) The numbers nw/ (x; — x;) by which one multiplies are generally of 
the order of magnitude of unity; hence it is easy to program their compu- 
tation for high speed machines. 

3) In computing u(x) from (1.8) the magnitudes of the first neglected 
divided differences in the region furnish an upper bound of the error, if | M| 
in (1.9) is no greater than unity (a case of frequent occurrence). 

4) The character of the differences is close to that of the familiar ordinary 
differences based on equal intervals between successive arguments. 


In employing modified divided differences over a large region X, it is 
often desirable to consider several subregions, in each of which w is chosen 
to be close to the average length of the interval over the subregion. 

Error Magnification. Consider the ordinary central difference of even 
order 2n, based on the entries u_», Un41, Uax—usually denoted by 
From (1.12) it is seen that the entry uw» appears in this difference with the 
coefficient (—1)"(%"). In the neighboring differences 6°", |k| <n, uo ap- 
pears with the coefficient (—1)"-*(2%,). It follows that if uo has an error €» 
which is substantially larger numerically than | 5?"u,| /(2"), and if the errors 
in neighboring entries are of a lower order of magnitude, then the differ- 
ences 5", will alternate® in sign, and their numerical ratios to &*u will 
be approximately (2",)/(%"). Since the error €9 is magnified most in 6"xo, 


_an incorrect entry can be picked out immediately, and even the order of 


magnitude of the error can be estimated from the difference pattern. (A 
divergence of the pattern from that of the binomial coefficient ratios indi- 
cates the presence of more than one error.) If a regular function u(x) is 
tabulated at equally spaced arguments, systematic differencing of the entries 
constitutes the most powerful tool for checking the function. 

The corresponding modified divided differences, based on unequal inter- 
vals between successive arguments, are inherently more difficult to interpret 
from the viewpoint of error indication. However, certain useful inequalities 
can be set down. In order to derive them, we shall limit the region X, and 
specify that the arguments x, form either an increasing or a decreasing 
sequence. In application this is actually the most important case; and in a 
region where x, — x,x_1 changes sign for some k, we may consider the sub- 
regions in which the sign of x, — x,-1 is constant. Furthermore, no essential 
restriction will be added if for convenience we consider only the case where 
X~ > X-1; this will now be assumed. 

Just as in the case of ordinary differences, u, enters into 2m + 1 con- 
secutive differences of order 2n. It will be enough to examine a fixed u» and 
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to study the coefficients with which it enters into the various differences of 
order 2n. Let us write 


(2.1) = [W; Xn] = MEP ur, 
k=—n 
and more generally 


(2.2) = (W; Xm = MSPs 
k=—n 


When no ambiguity is likely to arise, the superscripts will be dropped and 
we shall write for 
Let 


Wey = Xe — XE-1- 


Clearly c is positive, since x, — x,-1 is assumed positive. 
From (1.11) we have 


(2.3) Moo = (—1)"G")/Vao, 
where 
Co + + (co + + C2 at + 
Voo = Coc, 2 2 3 es 
x + Cx ( (2 +- 
Xe Xj 1 < =< 
(2.4) jo =1,ifk21. 
(x; — x0) c= -1ifk <0 
j=(n+le 
In particular 
(2.5) Mio — (xo — Xn) M_1.0 — (Xn — Xo) 


Moo (Xn41 — Moo (x9 — 
From (2.4) we have 


THEOREM 1. The coefficients My, of uo alternate in sign. 
THEOREM 2. If p< ce < P, for —n< k< n, then 


1 1 
n Dp — M10 n P 
(2.7) Moo ~n+1p’ 
2n 2n 
+ 


( < (= (, ») (fy. 


(2.8) () 


The above inequalities follow from (2.3), (2.4), and (2.5). 
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If P = p = 1 the inequalities give the well known relations among the 
coefficients of the ordinary central differences. From (2.6) and (2.8) we have 


THEOREM 3. If P < {2(2m + 1)/(m + 1)}}, the coefficients |M@| in- 
crease with n, for n = m. 


Proof. 
= | | w?(2m + 2)(2n + 1)_ | Mao|2(2n + 1) 
| (x0 — (Xn41 — 2 P?(n + 1) 
an +1 1 


since 2n + 1/(n + 1) is an increasing function of n. 

Coro.iary 3.1. If p> {2(2m + 1)/(m + 1)}4, | M&S | decreases with n. 

Before applying the previous results, it will be desirable to distinguish 
between three types of errors in the entries, just as in the case when ordi- 
nary differences are under consideration : 

Errors of Type (a). Systematic errors, prevalent in all entries over a 
region. Such errors cannot be discovered by differencing methods; hence 
they will be eliminated from the discussion. 

Errors of Type (b). Under this category we shall include 

1. An isolated error: An error in an entry—say in u»—which is con- 
siderably larger numerically than the errors of most entries in the region. 

2. Multiple errors: This term will be used when two or more entries in 
a region have errors which are considerably larger numerically than the 
errors of most entries in the region. 

The phrase considerably larger is necessarily qualitative, and depends 
on what degree of magnitudes the tests applied can distinguish. 

Errors of Type (c). 

1. Errors due to rounding approximate entries to a fixed number of 
decimals. 

2. Errors too small to be considered in Type (b). 

3. Errors of a non-systematic type which occur in too many entries of 
a region to be picked up by differencing tests. Usually all errors of Type (c) 
will occur in the last two decimal places of the entries. A discussion of errors 
of Type (c) will be given in Part II of this paper. 

Errors of Type (b). Consider an isolated error €o in uo. Let %o = uo + €0 
be the tabulated entry, and suppose that M§q’eo is of a higher order of 
magnitude numerically than 5° ux. 

Let 52'u, be the approximate value of the (2m)th divided difference in 
the region, and let 


(2.10) Soh” = — S = signum — 
On the assumption that Sg6” is due primarily to the error, we may write 


Sef”, 
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where p is a positive number close to unity. 


(2.11) hag | 
From (2.6) 
(2.12) < s 
The estimation of 52'u, for the region can be avoided. For we may study 
(2.13) deo — der = say. 
If €9 is an isolated error, then 
(2.13a) — Mio)eop = 
pS(—1)" => AAG 
From (2.7) 
< pS(—1)"o < (—1)"¢f? 
Maa(1 +5 7475): 
(2.14) 


THEOREM 4. If € is an isolated error, then 


Xo — Xn 
(2.15) — ~ — Xo 
to — by 1+ Xn — Xo 
Xo — 
Proof. 
2n M10 
— & €o(Moo — Mio) ). 

0, 0. 


— coMao( 1 - 
0,0 

Then (2.15) follows from (2.4) and the above. Since the right-hand side 
(2.15) is easy to evaluate, we have a criterion as to whether a possible error 
is an isolated one. The estimate of the error is furnished by (2.12) or (2.14). 

Let us now examine the effect of w on the location of errors. From the 
right-hand side of (2.10) it is clear that if w is chosen excessively small, 
then significant figures of u, will be moved off quickly to the right and it 
will be impossible to pick up significant errors of the last few places of the 
entries. Moreover, both P and # will tend to be large, since x, — x41 = qw 
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is fixed by the given grid for all k. In that case (2.14) will furnish a poor 
estimate of €9, due to the sharp variations in the inequalities. On the other 
hand, the choice of a very large w will simply have the effect of magnifying 
both small and large errors. A desirable w is one which tests as well as 
possible all the significant figures of the entries u%. For that purpose, w 
should be chosen as the average |x, — x:-;:| in the region. For then both 
P and p will tend to be of the order of magnitude of unity, provided suc- 
cessive intervals vary with reasonable continuity. 

If w is chosen as the average of the values x, — x,_1, then there may be 
some regions where P is relatively large with respect to w; and other regions 
where P is quite small. The question arises: If there is an error of Type (b), 
will (2.13) show it? From (2.13a) and (2.3) it is clear that gf will tend to 
be small in regions where the intervals x, — xs-1 are large with respect to 
the average interval of the space. Hence a small error will not show up. 
If the spread between P*" and ”* is small, in such a region then (2.14) can 
be used, to examine critically a fairly small g??. However, if an isolated 
group of k < m intervals is excessively large, then » and P will vary consid- 
erably, and the work of testing every doubtful entry may become too 
laborious. By the same considerations, a region where the values of x, — Xx_1 
are excessively small with respect to the average will show a big error mag- 
nification, and such errors will be readily traceable. Example 2 illustrates 
some features of error patterns. 

If P varies much from , the inequalities given here will not offer a good 
estimate of the magnitude of an isolated error. However, from Theorem 1 
it is assured that a sufficiently large error will show up, because of sign 
variation in the vicinity of an error, if the magnitude of the modified divided 
difference is of lower order than the error. Small errors, on the other hand, 
can be masked by the rounding errors in the vicinity, often to a greater 
degree than in the case of equally spaced entries. 

Let us now examine briefly the consequences of Theorem 3. In ordinary 
differences, where all c, are unity, the numerical value of Mo,o increases 
with m. Thus the magnification of an error increases in successive differences. 
If differences of high enough order are taken and the interval is sufficiently 
small for the correct differences to fall off numerically, an isolated error will 
be noticed more and more in successive differences. To have this property 
carry over to regions where all c, are not equal to unity, Theorem 3 requires 
that no c be numerically greater than 2[(m + 4)/(m + 1)]*. Thus in 
regions where the intervals between successive entries are numerically larger 
than twice the average interval in the region, successive differencing will not 
tend to magnify the error, and so detection of the error may become difficult. 
It would also be desirable to have the ratio | M+1,0/Mo,o| as small as possible. 
If the c are all equal to unity, this ratio is »/(m + 1) for the difference of 
order 2n. Hence with increasing 1, this ratio is so close to unity as not to 
offer much help in the variation of magnification even if all c, are unity. 
However, the fact that the sign of Mio/Moo ts always negative is of the 
utmost importance. For even if this ratio is numerically close to unity, the 
fact that successive M;,o differ in sign will help in detecting an error— 
especially so in regions where | Mo,o| increases with n. 

Although attention was fixed on even central differences, the main char- 
acteristics apply also to differences of odd order. 
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The following examples are instructive. 


Example 1 
(a) Unmodified Divided Differences of Yo(x) 

x Yo(x) 
1.00 .08825696 
1.04 11881884 __ 49959 
1.05 12624806 242922 _ “41953 S124 —.10 
1.08 114804063 2264190 —--1038 
1.20 22808350 670239 _ 086985 -0000 
1.25 25821685 8026670 _ "364579 -071895 05102-9540 
"5735327 05965 0390 
1:28 27542283 “3138327 _ "357021 95965 _'9307 = 02 
132 29736448 —"350673* —0326 
1.37 132321351 3169806 _"344473 04769 
1.41 34252663 4859780 -04343 
1.44 (35651952  -4022297 


Comment. In obtaining the unmodified divided differences, the computer 
was instructed to retain in [x,i41 --- Xi4n41] the same number of significant 
figures as there were in [xiXi41 --- Xn] — -- Xny1]. Another reason- 
able method is to retain in a column only the number of significant figures 
that are common to most entries in the column. The difficulty here is that 
the computer does not know, until the column is finished, what this common 


number should be. 
(b) Modified Divided Differences of Yo(x), in Units of the 8th Decimal Place 


w = 0.044 
x Vo(x) by? 6,5 

1.00 08825696 

1.04 11881884 «3361807 163509 
3268857 6338 

1.05 12624806 ~159751 
3196244 5305 198 

108 14804063 153318 ~720 
2934905 4446 119 

1.20 22808350 — 146582 —590 
2651735 3675 107 

125 (25821685 —141010 —459 —24 
2523544 3049 077 

128 [27542283 — 138239 —357 
2413581 2704 067 

132  .29736448 ~135781 —293 -9 
2274715 2437 059 

141 134265263 «3138303 134369 

144 "35651952 

Example 2 


Modified Divided Differences of x?, in Units of the Fifth Decimal Place 


w = 0.35 
—13.82400 
— 8.00000 —139106 
-15 — 3.27500 — 124950 
170100 
ay — 1.33100 — 71595 
78050 11434 
— [21600 — 44100 22231 
15050 40016 
‘1 — 00100 — 12250 —14291 
1050 20621 
2 ‘00800 14700 4203 
13650 25725 
5 "12500 = 93100 0 
‘6 "31600 — 200900 
7 "34300 117600 
1.0 1.00000 a 75950 
14 2.74400 


or 
or 
ig 
aS 
w 
th 
ns | 
), | 
to 
to | 
ip. : 
an 
ed 
id- 
00 
k—1 
ig- 
tes 
11 | 
ign x 
led 
nd, 
ter 
ses 
>es. 
tly 
will 
rty 
ires 
ger 
not 
ult. 
ble. 
e of 
t to 
lity. 
the 
the 
or— 
har- 


10 ON MODIFIED DIVIDED DIFFERENCES I 


Comments on Example 2. The entries corresponding to x = —1.5 and 
0.6 contain errors of the same magnitude. Yet the error in «(.6) is already 
evident from 6,?u, but not so the error in u( —1.5). This illustrates the fact 
that errors are magnified most in regions where the intervals c,w are small 
with respect to w. Let 


— x1 = — X% 
— Xo Xo — 


Xo — Xn Xn — Xo 
(1 + Xn+1 — =) /(1 + xo — 


From the table of differences, the left-hand side of (2.15) for this example is 
equal to 


x= .6=%, then 
using (2.15), 


—3.185 
—2.940 


Hence the error appears to be an isolated one. 

Let us now consider what information can be obtained from (2.12) about 
the nature of the error. In the usual case, the true value of 6@” is unknown, 
and gf” is estimated from adjacent differences. In the present case, the 
second differences at x = 0.1 and 0.2 warrant a guess that the second modi- 
fied difference at x = 0.6 is about 0.4. This is based on the assumption that 
the correct third modified difference in the region is about 0.25. Since, 
from (1.6), 


= 1.08 


[w; xoxrxe%s + [w; = [w; 


and (x3 — x9)/3w is approximately } in this region, we have, applying the 
above formula twice, 


6,2u = 0.15 + 0.25 = 9.4, at x = 0.6. 
Using the above estimate in (2.12) we arrive at the following calculations: 
P= .86; p=.285; P? = 0.74; p* = 0.081, 

| —2.009 — 0.4] = 2.409; & 1.205. 
Thus (2.12) implies 

1.205(.081) < |e] < 1.205(.74); or 0.097< |e] < 0.90. 
The true error does lie between the above limits, but the actual size of the 
error must be obtained from recomputation. The important fact to observe 
is that the pattern of the second differences correctly shows that (0.6) is 
in error. 


The error in u(—1.5) shows up in the fourth difference. 
Example 3. The following four values of f(y) are available. 


y f(y) A a 
3.7416573868 +.00824 2550 
3.777 +.00003 5971 — 230753 +79 
3.778 —.00019 4782 — 230674 
3.779 —.00042 5456 
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The differences of the equally-spaced values have been entered. It is 
required to find the zero of f(y). If it were certain that the three last values 
of f(y) are correct, and that A*y is negligible, this zero could be easily ob- 
tained by quadratic inverse interpolation. However, it happens that the 
computation of f(y) is quite laborious, and it is desirable to check the 
accuracy of f(y). Let us obtain modified divided differences, using the inter- 
val w = .001, since in this problem we merely want to have an indication 
of the magnitude of A?f in this region, at an interval of .001 in y. The values 
follow. 


w= .001 
y f(y) be 
3.7416573868 +.00824 2550 —232201 
3.777 +.00003 5971 330783 80 
3.778 —.00019 4782 —230674 79 
3.779 —.00042 5456 


The fact that the two values of 4,” differ by only one unit is assurance 
that the values of f(y) are correct, and that quadratic inverse interpolation 
is adequate. The solution is yo = 3.77715 586, to eight decimal places. [It 


e 
can be shown that A = }y¢° satisfies the system (+ ay t a) u = 0, 


u = 0 on the boundary C, u > 0 in interior of region where C is the ellipse 
+ = 1.) 


NBSINA 
Consolidated Engineering Corporation 
Pasadena, California 


lene L. M. MiLnE-Tuomson, The Calculus of Finite Differences. London and New York, 
2 The numbers u&, and uu: will be said to have the same order of magnitude if 
1/10 < |u:/u2! < 10. 
3 See J. C. P. Mrtxer, “Checking by Differences I,” MTAC, v. 4, 1950, p. 3-11. 


GERTRUDE BLANCH 


A Minimum Problem Solved by Mesh Methods 


Introduction. In the following, the function y(x), subject to y(0) = y(1) 
= 1, is sought which will minimize the integral 


(1) f + de. 


This problem can be solved by the usual methods of the calculus of 
variations! but the differential equation involved is rather complicated. 
It is proposed here to solve the problem by mesh or “‘assumed polynomial” 
methods.” Methods of this sort have been rather extensively applied to the 
solution of differential equations and boundary value problems, and, re- 
cently, also to the determination of characteristic numbers or eigenvalues.* 

Basic Equations. To solve the above problem by mesh methods it is 
first assumed that the interval 0 < x < 1 is divided into segments each of 
equal length h = 1/(2m) by the points x;, ¢ = 0(1)2mn. It should be noted 
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that the solution will be symmetric about x = 4; hence only the portion 
0< x< 3} need be considered together with the condition 


(2) = 0. 


We then write, using approximate integration, 


(3) I = 2% 2hri/ (i + 
i=] 
where 
(4) r= {1+ — 
so that J is now a function of y; for 7 = 0(1)”. At a minimum we must have 
ol 
(5) 
If we write 
Se = Vi + Ve-1, Be = + 


(Sere), 


we can write equation (5) as 


Pk 


(6) = Ve + pe (Ye — Ye-1) + 


The process of solution followed here is to assume a curve y(x) and use 
this to determine the p’s and e’s at mesh points. Then (6) provides a recur- 
rence relation for determining new values of yo, ys, ---, ¥n which we write 
as Yo, Y3, ---, Y, since the mid-point slope condition (2) is not as yet satis- 
fied. Using second order deriving coefficients as given by SOUTHWELL’ or 
Mine! (p. 96-98) we calculate Y,’ which usually is not zero. Then, letting 
yi = Y/ — Y,’ so as to satisfy y,’ = 0 we integrate y,’ by use of inte- 
grating coefficients? to get a new set of y’s. The above process is then 
repeated. 

The accompanying table shows results for three trials using h = 1/8 
and starting with y = 1. An approximate value of 3] was found for each 
assumed curve y(x). These successive values, i.e. .5, .48148, .48144, indi- 
cate a minimum being reached. 


Yk Tk Sk Pk 8k ek Y;’ Yk 
1.00000 — 1.00000 +.0625 1.00000 
1.00000 .00000 1.00000 2.00000 1.00000 — — 1.00000 —.0625 1.05469 
1.00000 .00000 1.00000 2.00000 1.00000 .50000 —.015625 .984375 —.1875 1.09375 
1.00000 .00000 1.00000 2.00000 1.00000 .50000 —.015625 .953125 —.3125 1.11719 
1.00000 .00000 1.00000 2.00000 1.00000 .50000 —.015625 .906250 —.4375 1.12500 


k 
0 
1 
2 
3 
4 
0 1.00000 — 1.00000 .50504 1.00000 
1 1.05469 .05469 1.09152 2.05469 1.003638 — _ 1.05469 .37000 1.05362 
2 1.09375 .03906 1.04769 2.14844 .99937 .48553 —.01708 1.09250 .24124 1.09071 
3 1.11719 .02344 1.01743 2.21094 .99870 .43512 —.01529 1.11500 .12268 1.11241 
4 1.12500 .00781 1.00195 2.24219 — 40744 —.01430 1.12347 .00804 1.11953 
0 
1 
2 
3 
4 


1.00000  — 1.00000 .49656 1.00000 
1.05362 .05362 1.08812 2.05362 1.00096 — 1.05362  .36136 1.05358 
1.09071 .03709 1.04309 2.14433 .99971 .48486 —.01695 1.09034  .23280 1.09064 
1.11241 .02170 1.01496 2.20312 .99976 .43596 —.01523 1.11182 .11448 1.11234 
1.11953 1.11947 


.00712 1.00162 2.23194 —  .41018 —.01433 1.11896 —.00024 
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Conclusions. It can be seen from the table that the process converges 
with fair rapidity. Further calculations by the author indicate that the 
results obtained after three trials contain an error of less than .001 for 
y(}). If the usual methods of the calculus of variations were employed the 
resulting non-linear differential equation would presumably have to be 
solved by finite difference methods anyway and it does not appear that this 
would be as easy a computation to carry through as the above. 

In writing (3), first order divided differences have been used, these 
being the simplest and at the same time adequate. Higher order expressions 
for the derivatives may be employed but will in general result in more 
complicated recurrence relations. The iterative procedure for the solution 
apparently must be devised anew for each different class of problems. 


G. C. BEst 
3717 Crestline 


Fort Worth, Texas 
1R. Wernstock, Calculus of Variations. New York, 1952. 
2R. V. SouTHWELL, Relaxation Methods in Theoretical Physics. Oxford, 1946. 


Cottatz, Numerisches Behandlung von Differentialglerchungen. Berlin, 1951. 
4W. E. Mitne, Numerical Calculus. Princeton, 1949. 


A Modification of the Aitken-Neville Linear 
Iterative Procedures for Polynomial 
Interpolation 


A. C. AITKEN" has described a method of interpolation which is equiva- 
lent to the use of Lagrange’s polynomial formula but consists principally 
of the repeated computation of an attractively simple algorithm, well suited 
to desk calculators. The method does not require uniform spacing in the 
values of the argument at the points at which the values of the required 
function are given, though uniformity permits of a convenient check of some 
aspects of the calculation. It can therefore be used for both direct and 
inverse interpolation, and it is particularly valuable for the latter. The 
procedure depends on the following property, which is also discussed at 
greater length by E. H. NEvILie? and by W. E. Mrtne? 


On the basis of the known values of a function u at m values of the 
argument X—that is, at a point P and m — 1 other points, denoted col- 
lectively by Q—we may have obtained a polynomial interpolate upg of de- 
gree n — 1 for the value of the function at some other point X = x. We may 
have obtained also another interpolate ugar of degree n — 1, for the same 
value x, on the basis of the » — 1 points Q and a further point R. Then the 
polynomial interpolate of degree n for X = x, based on all m + 1 points P, 
Q, R, is 


In this formula, xp and xz are the values of the argument at the points 
P and R respectively. It is easily computed on most desk calculators, espe- 


> 
e 
g 
n 
8 

00 | 
69 
75 
19 
00 
00 
62 
71 
41 
153 

upg Xp 
u = Xr — Xp). 

58 
164 
134 
47 
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cially as the divisor xz — xp can be found (or checked), without extra work, 
as the net total of the two multipliers xz — x and —(xp — x) which are 
used in evaluating the determinant. Moreover, the linearity of the formula 
permits upg and ugp to be adjusted equally by any suitable amount before 
interpolation, thus minimizing the number of digits to be used at each stage. 
We use the term ‘interpolation’ loosely here to include extrapolation, which 
strictly occurs when xp — x and xz — x have the same sign. Aitken! calls 
this elementary process ‘linear iteration by proportional parts,’ and has 
also (according to WoMERSLEY‘) called the derived value of upgz a ‘linear 
cross-mean’ of upg and ugr. 

This iterative linear process or algorithm can therefore be used to in- 
crease the order of approximation of an interpolated value for a function, 
starting with a simple linear interpolation between two known values 
(m = 1) and progressively incorporating others. The final result, when the 
same known values have all been used, will necessarily be the same (apart 
from rounding-off variations) regardless of the order in which they are 
incorporated. The algorithm will be computed 4n(n — 1) times to incor- 
porate m points, but with a suitable schedule considerable economies and 
simplifications can be effected. Aitken! proposes a schedule in which the 
work sheet for the computations can be laid out in such a way that the 
values required for each determinant always appear at the corners of rec- 
tangles. Neville? proposes a schedule which lacks this advantage but which 
usually gives a more rapidly convergent process. Neville’s schedule permits 
of the addition of further known points, if desired, after the completion of 
the calculation of an interpolate, with greater freedom than Aitken’s does, 
owing to the fundamental asymmetry of the latter. As Neville? and K1n- 
carp‘ have demonstrated, Neville’s method is also readily adapted to utilize 
the known values of differential coefficients at points at which the function 
is known. 

The convenience of Aitken’s schedule as regards the determinants, and 
much of the convenience of Neville’s as regards the latitude in the addition 
of further points, are possessed by a third schedule which is also usually a 
little more rapidly convergent than Neville’s in the early stages. With a 
condensed notation, using (7, j, k, ---) to stand for the interpolate based on 
the points X = x;, x;, X,, «++, the new procedure follows the schedule shown 
diagrammatically below for the case m = 6. The schedule for smaller values 
of is the appropriate fragment of this one, and the procedure generalizes 
immediately to greater m. The cross products involving each line are taken 
with the items nearest the axis of symmetry of the schedule, and alternately 
on the further and nearer sides of this axis for successive columns, as is 
indicated by the lines in the diagram. In this diagram, the value of each 
interpolate (7, j,k, ---) is written on the same line as one of the values on 
which it is directly based. Certain entries, marked by asterisks, are dupli- 
cated so that the values used in the determinant always occur at the corners 
of a rectangle. These entries can be computed separately as a partial check 
without interrupting the rhythm of the process. Incidentally, the geomet- 
rical interpretation given by Aitken! (p. 73-74) for his schedule can be 
adapted to represent the new schedule, and also to represent Neville’s. 
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(1 1,3 1,2,3,4,5 1,2,3,4,5,6) <x 
(2) 2,3,4 2,3,4, x 
(3 
(4 
(S 4,5 2,3,4,5)8 
(6 3,4,6 1,2,3,4,5,6)* 


The new schedule can be regarded as a modification of Aitken’s to 
introduce properties of symmetry which are often—though not always— 
desirable in an interpolation formula. Both for ease of computation and for 
accuracy of interpolation of a given order, it is normally most satisfactory 
when the numbers of given points on each side of the required point are 
either equal or differ by unity. That is to say, the required point should lie 
in the interval enclosed by the two given points which lie nearest to the axis 
of symmetry of the general pattern. Then the linear cross-mean procedure 
used at each stage is equivalent to an interpolation (for the columns of 
interpolates of odd degree) or an extrapolation (for the columns of inter- 
polates of even degree) on a straight line fitted as a chord to a curve repre- 
senting the previous column of values, and this concept is useful in judging 
the closeness of approximation in the new method. Moreover, the latter is 
designed to make the fullest use of the points nearest the required value 
from the earliest stages onwards, and to use them for interpolation rather 
than for extrapolation. Both of these features normally increase the rate of 
convergence, and thus facilitate both the computations and the detection 
of errors. 

In the six-point example used by Milne® to illustrate Aitken’s and 
Neville’s methods, the range of variation of the results in the earlier columns 
is considerably reduced under the new schedule. We shall not, however, 
use this as an example here, since the given values of the function are at 
uniformly spaced values of the argument, and the NBS Tables of Lagrangian 
Interpolation Coefficients* enable the final interpolate to be found in two 
minutes easily with a calculator. The time taken to compute the fifteen 
linear cross-means under Aitken’s or Neville’s or the new schedule cannot 
practicably be reduced much below ten minutes (including writing down 
the necessary digits of the intermediate interpolates) with a fast electrical 
desk calculator. 

WHITTAKER & RoBInson’ demonstrated a method of inverse interpola- 
tion for finding the positive root of the equation u’? + 28u‘ — 480 = 0, 
given the values of x = u? + 28u* — 480 at the five points wu = 1.90(.01)1.94. 
Aitken! (p. 71) and Neville? (p. 94-97) used this example to illustrate their 
schedules of linear iteration by proportional parts. Working to 10D for u, 
with 8D for x, the new schedule may be applied to solve this problem as 
follows: 


u=1.9 22 88 41 

1.90 2271 69289 88 48558 41731 527 —25.71402 610 

1.91 2279 04030 88 45462 41643 —14.62541 674 

1.92 2286 32562 — 3.30746 392 

1.93 2286 32562* + 8.24394 354 

1.94 2283 41050 88 36413 41643* +20.03298 301 
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Thus the estimated root is u = 1.92288 41527, which differs by six units 
in the last place from the correct 10D value, of which the last two digits 
are 33. The complete inverse interpolation process takes about 12 to 15 
minutes on the desk calculator referred to above. Further values of x;, say 
for u; = 1.95 first, followed by u; = 1.89, 1.96, 1.88, 1.97, etc., could be 
added if desired, and the existing results incorporated completely into the 
extended schedule up to the point at which the number of decimal places 
carried from the start ceases to justify it. In this example not more than 
two values can advantageously be added. 


M. C. K. TWEEDIE 
Statistical Laboratory 
Virginia Agricultural Experiment Station 
Virginia Polytechnic Institute 
Blacksburg, Virginia 


1A. C. ArTKEn, “On interpolation by iteration of proportional pa without the use 
of oS a Edinbur h Math. Soc., Proc., ser. 2, v. 3, 1932, p. 56-76 
. 2E. H. NEvILte, “Iterative interpolation. ” Indian Math. Soc., Jn., v. 20, 1933, p. 
87-120. 


3W. E. Mine, Numerical Calculus. Princeton University Press, 1949, chap. III. 
0-117 R. WoMERSLEY, “Scientific computing in Great Britain.” MTAC, v. 2, 1946, p. 
110-117. 
5 W. M. Kincanp, “Solution of equations by interpolation.” Annals of Math. Statistics, 
v. 19, 1948, p. 207-219. 
°'N BSMTP, Tables of Lagrangian Interpolation Coefficients. New York, Columbia Uni- 
versi Press, 1 1944. 
*, = WHITTAKER & G. Rosinson, The Calculus of Observations. Fourth ed., London, 
1948, chap. IV. 


RECENT MATHEMATICAL TABLES 
1144[F ].—J. Toucuarp, ‘‘On prime numbers and perfect numbers,” Scripta 


Math., v. 19, 1953, p. 35-39. 
This note contains a small table of the function 


C,(n) = k?a(k)o(n — k) 


for p = 0(1)3 and m = 2(1)11. Here o(k) denotes the sum of the divisors 
of k. These coefficients occur in the expansion of certain elliptic functions. 
More specifically, 


X Cy(m)x™ = pPo1 
n=1 
where 


&.= > m'n*x™. 


n,m=1 


D. H. L. 


1145[1].—E. L. Kaptan, ‘Numerical integration near a singularity,”’ Jn. 
Math. Phys., v. 31, 1952, p. 1-28. 


The types of singularities considered in this paper are of the forms 


x4A(x), xtA(x), x"{A(x) logx + B(x)], n =0,1, 
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where A and B are regular near the origin. There are 14 principal tables of 
coefficients ranging from 3 to 6-term formulas. The accuracy is from 10 to 
14D. Certain auxiliary tables of elementary symmetric functions of sets 
of 2, 3, or 4 small integers are included. The paper suffers from lack of 
illustrative material. The same problem for x+A (x) has been considered by 
Y. L. Luxe [MTAC, v. 6, p. 215-219], who discusses also the paper under 
review. 


D. H. L. 


1146[I].—Taxaniko YAMANOUCHI, Tables of Coefficients of Everett's Inter- 
polation Formula. Report no. 1, Computation Institute, Tokyo, 1946, 
19 p. 15 X 21 cm. 


These tables give coefficients A and B for Everett's interpolation formula 
y (xo + Oh) = Oy: + — (Ady, + Bé*y0) 
A=}0(1-@) and B= -¢), o=1-8, 
and for the inverse interpolation formula 
6, = n + aA (0,1) + 6B (@,-1), 


where 


where 
n = 4% (y — — Yo), = — Yo), 
B = — and x = + O,h. 


The coefficients A and B are tabulated for @ = 0(.001)1 [alsog = 1( —.001)0] 
to 5D, with differences and proportional parts. 


H. E. SALzer 
Diamond Ordnance Fuze Laboratory 
Washington 25, D. C. 


1147[K].—R. A. BrapLey & M. E. Terry, “Rank analysis of incomplete 
block designs. I. The method of paired comparisons,” Biometrika, v. 39, 
1952, p. 324-345. 


Given ¢-treatments with true non-negative preference ratings, +; and 
estimates, ~:; (7: = Sp; = 1). The treatments are compared in blocks 
of two, (3) blocks constituting a basic design; m repetitions of the basic 
design are considered. A likelihood test statistic (B) is given to test the null 
hypothesis, Ho: 7; = 1/t, against alternative hypotheses in which m groups 
of treatments can be formed so that the treatments differ from group to 
group but are alike within groups. Two special cases are considered: (i) 
m = t, (ii) m = 2. 

Tables (Appendix A) are prepared to make significance tests for case (i), 
when ¢ = 3 and 4. For ¢ = 3, m = 1(1)10 and for ¢ = 4, m = 1(1)6. These 
tables are in terms of all possible combinations of rank sums for each treat- 
ment, where 7;;, (= 1 or 2) is the rank of the 7-th treatment in the k-th 
repetition of the block in which both treatments i and j appear. Values of 
the following are given for each combination of rank sums: p; to 2D for 
4 = 1(1)t; the test statistic, B, to 3D; the significance probability, P (proba- 
bility of this or a smaller value of B,) to 4D. 

A combined analysis is also considered for g groups of rank sums in which 
the x’s may change from group to group, so that it is not desirable to pool 
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the rank sums. The test statistic B,° is found by adding values of B, for the 
individual groups. Appendix B gives values of B,* to 4D and P to 3D for 
n = n'g,n', g = 2(1)5. 
ta The authors also set up a test of agreement of the rank sums from group 
to group, and large sample tests for this agreement and for B,, Bz and B;,°. 

An example is presented for ¢ = 3, g = 2, n’ = 5. 

R. L. ANDERSON 

North Carolina State College 
Raleigh, North Carolina 


1148[K ].—H. E. Dantets, ‘‘The covering circle of a sample from a circular 
normal distribution,” Biometrika, v. 39, 1952, p. 137-143. 


If a sample of observations is taken from the circular normal dis- 


tribution 
dF = exp [ — (x? + dxdy 


the covering circle is defined as the least circle (radius r and centered at 
distance p from the mean) such that all the observations are on or within 
the circle. The joint distribution of r and p is given by 


dF, (r, p) = —n(n — 1)o*{exp [—(r? + p*)o* p) drdp 


where 


P(r, 6) = exp ds 
and 
I,{s) = (24 =f [exp (z cos ] cos wi dé. 


The distribution for r taken singly is found to be 
F,(r) = n{1 — exp — (m — 1){1 — exp 


A table is exhibited showing the mean to 3D, the variance to 4D, and the 
1%, 5%, 95% and 99% points to 3D of the distribution of r/c for samples 
of m = 2(1)15, 20(10)50, 100. 

T. A. BICKERSTAFF 
University of Mississippi 
University, Mississippi 


1149[K].—H. A. Davin, “Upper 5 and 1% points of the maximum 
F-ratio,”’ Biometrika, v. 39, 1952, p. 422-424. 


The ratio of the largest to the smallest of k sample variances, all with 
the same number of degrees of freedom », was introduced by Hartley' as a 
test statistic for the hypothesis that the population variances are equal, 
under the assumption of normality. Hartley gave a table of the upper 5% 
points for this ratio, containing values which were exact for k = 2 but were 
obtained from an approximate expression for k > 2. In the present paper 
numerical quadrature is used to evaluate the exact upper 5% and 1% points 
for this ratio with an accuracy which may possibly leave the third digit in 
doubt. Tables thus computed are presented for k = 2(1)12, » = 2(1)10, 12, 
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15, 20, 30, 60, ~. It is found that Hartley’s approximate table tended to 
underestimate the 5% points. 

Z. W. BIRNBAUM 
University of Washington 
Seattle, Washington 


1 iH. O. Hart ey, ‘The maximum F-ratio as a short-cut test for heterogeneity of vari- 
ance,”’ Biometrika, v. 37, 1950, p. 308-312 [MTAC, v. 5, p. 145]. 


1150(K ].—M. H. Gorpon, E. H. Lovetanp, & E. E. Cureton, “An ex- 
tended table of chi-square for two degrees of freedom, for use in com- 
bining probabilities from independent samples,” Psychometrika, v. 17, 
1952, p. 311-316. 


These tables give the value of chi-square for two degrees of freedom 
corresponding to a given probability value. The values of chi-square are 
given to 4D for the argument, » = 0(.001).999. The major function of 
these tables, according to the authors, is to convert probabilities obtained 
from independent samples into corresponding values of chi-square so that 
they may be combined by the technique proposed by Fisher if one is inter- 
ested in making a new test of significance of the combined data. A numerical 
illustration showing the application of the Fisher process and the use of the 
tables accompanies the article. 

C. F. Kossack 
Purdue University 
Lafayette, Indiana 


1151[K ].—E. J. GumBEL, J. ARTHUR GREENWOOD, & D. Duranp, “The cir- 
cular normal distribution: theory and tables,’’ Amer. Stat. Assn., Jn., 
v. 48, 1953, p. 131-152. 


Let a1, a2, -++, a, denote the angles determined by points on the cir- 
cumference of a unit circle, with respect to a fixed radius, and let a» denote 
the angle determined by their centroid, whose radius vector is @. Then ap 

n 


is a species of mean given by > sin (a; — ao) = 0. In this paper the proba- 


i=1 
bility density function f(a), for which the maximum likelihood estimate 
of its location parameter is a mean of the form ap, is found to be 
f(a) = {exp [k cos (a — ap) ]}/[2aJ0(k) ], where J is the Bessel function of 
the first kind of purely imaginary argument. This is analogous to the ordinary 
normal distribution in the case where the ” points are situated along a straight 
line instead of a circle. Table 2 gives the maximum likelihood estimate of 
the scale parameter k determined from the relation Jo’(k) — @Jo(k) = 0, 
for the radius vector & for the given observations. Values of k are given to 
5D for @ = .00(.01).87 with second central differences that are modified 
when @> .69. 

Table 3 gives 3D values for the function ¥(a) = {[exp (k cos a) ]/Io(k)}! 
for k = 0(.1)4 and a = 0°(10°)180°. This function is used in equi-areal (i.e., 
area-preserving) plotting of the circular normal distribution f(a). 


Table 4 gives 5D values of ®(a) = f ; f(x)dx for a = 5°(5°)180°, 
k = 0(.2)4, together with the second central differences in each direction, 
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6.7, 5,2. This table may be used for finding areas or significance levels. It is 
an abridgment of a 7D table, with intervals half as large, that was computed 
from the Fourier cosine series for exp (& cos x), and from tables of I,,(%) to 
n = 14 obtained (by arrangement) from the proof sheets of Part II of the 
BAASMTC’s Table of Bessel Functions.' 
J. LiEBLEIN 

National Bureau of Standards 
Washington, D. C. 

1BAASMTC, Bessel Functions, Part II. Cambridge, 1952 [MTAC, v. 7, p. 97]. 


1152[K ].—A. K. Gupta, “Estimation of the mean and standard deviation 
of a normal population from a censored sample,’ Biometrika, v. 39, 
1952, p. 260-273. 


The author recognizes two ways of censoring a sample: (1) observations 
above or below a given (truncation) point may be censored (type I); and 
(2) the n — k smallest or greatest observations may be censored (type II). 
He derives maximum-likelihood equations for estimating the mean and 
standard deviation of a normal population from a type II censored sample. 
Although differing in algebraic form, the author’s results when applied to 
any given sample with the largest (or smallest) measured sample observa- 
tion, x,, taken as the truncation point, lead to identical estimates as corre- 
sponding estimating equations derived earlier by STEVENS,’ HALD,? and the 
reviewer.* The author indicates that samples considered in the above refer- 
ences are of type I since there, the truncation point is assumed known prior 
to collecting sample data. It must, however, be pointed out that estimating 
equations in each of these papers were derived with the total number of 
sample observations, m, the number of measured observations, k, and hence 
the number of censored observations, 2 — k, considered as known constants. 

Tables are provided in this paper for computing estimates when the 
right (larger) tail of the sample is censored. The function, 2(y, »), used in 
the estimation process, is tabulated to 4D for y = .05(.05).95 and for 


_p = .1(.1)1.0, where 


= s*/(s? +d’) = (1 + nz — 2)/(1 + 92) 


and 
fexp (—1°/2)} exp (—f/2) a) 


in which s? is the sample variance, d = x, — Z, x; is the largest measured 
observation, Z is the sample mean, p = k/n, and » = (xx — u)/o, where wu 
and o@ are the population mean and standard deviation. For computing 
asymptotic variances and covariances of the estimates, tables of o;; (¢, 7 = 
1,2) are given to 5D for p = .05(.05).95(.01).99, where [o.;] = [viz]“ 
in which »,; are elements of the maximum likelihood variance-covariance 
matrix. 

Best (minimum variance) linear estimates, u* and o* from type II 
samples for small m are also considered in this paper, and a table with 
entries to 5D of the coefficients 8; is provided, where 


= b= 2(1)e —1; = 3(1)10 
i=] 
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in which the are the x;,4 = 1, ---, k, rearranged 80 <Xejn- 
A similar table with entries also to 5D is given of the coefficients y;, where 


k 
o* = 2(1)n; = 2(1)10. 
i=l 


Two final tables with entries to 5D are given of the variances of estimates 
o* and u* for ,k = 2(1)10. 
A. C. CoHEN, Jr. 
University of Georgia 
Athens, Georgia 
1W. L. Stevens, “The truncated normal distribution” (Appendix to paper by C. I 


aes — calculation of the time mortality curve’), Ann. Appl. Biol., v. 24, 1937, p. 


2 A. HALp, ‘‘Maximum likelihood estimation of the parameters of a normal distribution 


which is truncated at a known point,” Skandinavisk Aktuarietidskrift, v. 32, 1949, p. 119- 
134. MTAC, v. 5, p. 74 


3A. C. CoHEN, jr., “Estimating the mean and variance of normal populations from 
singly truncated and doubly truncated samples,” Ann. Math. Stat., v. 21, 1950, p. 557-569. 


1153[K ].—M. R. Samprorp, ‘‘The estimation of response-time distribu- 


tions. II]. Multi-stimulus distributions,’ Biometrics, v. 8, 1952, p. 307- 
369. 


“The types of response time distribution occurring when two or more 
stimuli act on a sample of individuals are discussed. Two situations are 
discussed in considerable detail ; the ‘accidental death’ model, in which each 
response can be related to its appropriate stimulus, and the ‘natural death’ 
model in which the exact cause of any death cannot be determined, but 
the distribution of ‘potential survival times’ to the two stimuli can be 
assumed bivariate normal. Maximum likelihood methods of estimation are 
developed for these situations, and tables are given to simplify the cal- 
culations.” (From the author’s summary.) Tables of », A, and & are given 
to 5D for = —5(.1) — 3(.01)3(.1)5, where Z(n) = (2x)? exp (—7°/2), 


should be noted that é is always positive. In the tables a curious misprint 
occurs for all 7 = —2.75 in attributing to £ a minus sign. 


L. A. AROIAN 
Research and Development Laboratories 


Hughes Aircraft Co. 
Culver City, California 


1154[K ].—M. E. Terry, “Some rank order tests which are most powerful 
against specific parametric alternatives,’ Annals Math. Stat., v. 23, 
1952, p. 346-366. 


A statistic ¢,(R) based on a method of HoEFrrp1nc! yields a most power- 
ful rank order test of Ho: two samples of m and m observations come from 
the same continuous population, against the alternative H,: the two samples 
come from two normal populations with common variance o*, and with 
means 6 and ¢, respectively, where (6 — ¢)/o is positive and small. An 
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approximation to the null distribution of c; in terms of the ¢-distribution 
and a comparison of the exact and approximate critical values for N = 
6(1)10 for certain m and n is given. Under certain conditions, the asymptotic 
distribution of c,; under Hy is shown to be normal. Table 1 gives the exact 
null distribution of ¢,(R) to 2D for N = 2(1)10 for m< n, m+n=WN 
and distinct permutations of R (i.e. permutations of m 0’s and m 1’s), to- 
gether with the corresponding values of the MANN-WHITNEY’ U test. 
Table 2 gives the probability p of exceeding the critical values c; under Ho 
for all p< .1 to 3D for N = 6(1)10. 
INGRAM OLKIN 

Michigan State College 
East Lansing, Michigan 

1W. HoeErrp1ne, “ ‘Optimum’ nonparametric tests,” Proceedings of the Second Berkeley 
= pr —— Statistics and Probability, Univ. of California Press, Berkeley, 


2H. B. Mann & D. R. Wurtney, “On a test of whether one of two random variables is 
stochastically larger than the other,” Annals of Math. Stat., v. 18, 1947, p. 50-60. 


1155[K ].—J. E. Watsu, ‘‘Operating characteristics for tests of the stability 
of a normal population,’’ Am. Stat. Assn., Jn., v. 47, 1952, p. 191-202. 


Let m independent observations be drawn from a normal population 
with unknown mean m’ and unknown standard deviation o’. The problem 
is to test whether a small sample came from a specified normal population 
with known mean m and known standard deviation c. In this situation, the 
null hypothesis tested is that m’ = m and o’ = co. 

The most useful result of the paper is a set of tables which give the 
operating characteristic (OC) function values for the following three types 
of tests. The first test, and the most common, is based on the “limits” for 
the sample mean, , and the standard deviation, s (using m — 1). In the 
second test, let be replaced by the #-statistic 


t = ni(z — m)/s, 


so that the limits for ¢ and s are used. In the third test, the statistic s is 
replaced by 


s(m) = (x; 


to give a test based on the limits of and s(m). 

Table 1 provides ‘‘limits,” given to 3D, which define the critical region 
for ¢ and for s at a = .005, .01(.01).05 and m = 3,5. Table 2 provides 
“limits,”’ given to 3D, which define the critical region for and for s(m) at 
the same values of a and 7 as given in Table 1. If one or both of the statistics 
fall outside these limits, the null hypothesis is rejected. 

The OC function for a test is defined to be the probability that both 
statistics have values within their limits given the true values of m’ and o’. 
However, in this paper, the OC function values are given for alternative 
hypotheses which are expressed in terms of two equivalent parameters @ 
and b, where 


a =ni(m —m')/c’, b = a/c’. 
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Table 3 contains an OC function comparison, given to 3D, for the three 
types of tests at values of a = .01;” = 3,5;a = 0(1)4; 5 = 1/8, 1/4, 1/2, 
1, 2, 4, 8. Table 4 contains the OC function values, given to 3D, for the test 
based on Z and s at values of a = .001, .005, .01, .02, .05; m = 3, 5, 7; 
a = 0(1)4; b = 1/8, 1/4, 1/2, 1, 2, 4, 8. 

The results of this paper are directly applicable to those quality control 
situations for which the assumptions underlying the test are satisfied. This 
reviewer would enjoy the extension of the OC function values to the and 
R test together with the three types of tests included in this paper. 


G. W. McELratu 
University of Minnesota 
Minneapolis, Minnesota 


1156[L].—S. CHANDRASEKHAR & Donna ELsert, “The roots of 


(Am) — (A) = 0,” Cambridge Phil. Soc., 
Proc., v. 49, 1953, p. 446-448. 


The authors tabulate to 6S the first root, \;, of the equation mentioned 
in the title, together with 6S or 7D values of 


2 FR 


for » = .2,1 = 1(1)5; 9 = .3,/ = 1(1)6; and » = .4, .6, .8, 7 = 1(1)15. 
A. E. 


1157[L].—W. J. Duncan, Normalised Orthogonal Deflexion Functions for 
Beams. Aeronautical Research Council Reports and Memoranda, No. 
2281. His Majesty’s Stationery Office, London, 1950. 23 p. 


The orthogonal deflexion functions are 


4 
Aa) = lena — psy, 


and can be shown to be numerical multiples of (1 — x*)*P,,’’"(x) where 


x = 1 — 2&, m = 2n for S,, m = 2n + 1 for Ax, and P,, is the Legendre 
polynomial. 


The Green’s functions for a doubly built-in uniform beam are 


= Gilt, 8), = 2 GG, 9). 


Tables 1 to 5 were computed under the direction of R. A. FAIRTHORNE, 
and give 5D values of S;, S2, S3, A1, As for — = 0(.01)1. 

Tables 6 and 7 were computed by S. Kirxsy and give 5D values of 
G, and G, for &, r = 0(.05)1. 


A. E. 
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1158[L ].—Costantino Fasso, “Di un integrale intervenuto in una ques- 
tione di idraulica,’”’ Ist. Lombardo Sci. Lett., Rend., Cl. Mat. Nat., s. 3, 
v. 15, 1951, p. 471-497. 


The author denotes by ®(x, z) a suitably defined indefinite integral of 


= 
z+x 


with respect to x. He gives 5D values of ®(x,z) for x = .86, .9(.02)1.1, 
1.14 and z = 10, 5, 2.5, 1.5, 1, .5, .25, .1, 0 and for fewer x values and 
—z = .25, .5, .7, 1, 1.2, 1.5, 2.5, 5, 10. He also gives some numerical data 
and diagrams of the level curves of the surface, 


in (x, y, 2) space. 
A. E. 


1159[L].—Marion C. Gray, “Legendre functions of fractional order,” 
Quart. Appl. Math., v. 11, 1953, p. 311-318. 


The author gives a 6D table of P,(cos @) for v = .1(.1)2, @ = 10°(10°) 
170°, 175°. In comparison with the tables described in RMT 1110, the 
present table is much shorter in the » direction, much thinner in the @ 
direction, and there are no associated Legendre functions; but on the other 
hand there is one more decimal, and the tabulation extends to @ = 175° 
(instead of @ = 90°). 

There is a discussion of the zeros of P,(cos @) for fixed @ and variable v 
(see also RMT 1120 and the literature quoted there), and there are several 
useful expansions for numerical computation. 

A. E. 


1160[L, V].—J. BARKLEY Rosser & R. J. WALKER, Properties and Tables 
of Generalized Rocket Functions for Use in the Theory of Rockets with a 
Constant Slow Spin. Cornell University, 1953, 21.5 X 28 cm. 


This report treats certain indefinite integrals that describe the motion 
of fin-stabilized rockets having a constant slow spin about the axis of sym- 
metry. The integrands are products of trigonometric functions and Fresnel 
integrals. In terms of functions defined in MTAC, vol. 2, p. 213 and MTAC, 
v. 3, p. 474, the integrals tabulated are 


Rrc(a, x) = ft Rr(w) cos a(w — x) dw 


and the similar integrals obtained when Rr(w) is replaced by Ri(w) and 
when cos is replaced by sin. About half the space is devoted to a tabulation 
of these four integrals to 5D for x = 0.0(0.1)5.0 and a/w = 0.1(0.1)2.0. 
For a lesser range of x, the tables extend over a/x = 2.0(0.1)8.0. Second 
differences with ‘‘throw-back”’ are given for x. Auxiliary formulas and tables 
enable the entire range of real a and x to be covered. Careful attention is 
given to the accuracy attainable in the various regions. About half the 
space is devoted to definitions of the various rocket functions, to their 
connection with differential equations, to a review of the properties of the 
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simpler Rr and Ri functions referred to above, to a development of the 
properties of the functions tabulated, and to the methods of computation 
of the tables. 

LEvERETT Davis, JR. 
California Institute of Technology 
Pasadena, California 


1161[L ].—Tatupiro Sasakl, Tables of u(x) = Rp (ix) and oo(x) =arg T'(1—ix). 
Numerical Computation Bureau, Tokyo, 1950. 9 pages, mimeographed. 
These tables contain values of u(x) = Rey (ix) and oo(x) = arg I'(1 — ix) 

where x is a positive variable, ['(z) is the gamma function and ¥(z) is its 

logarithmic derivative. The values of u(x) are tabulated to 6D with second 
central differences for x = 0(.01).5(.02)2(.05)2.5. The table of oo(x) is also 
to six decimal places with second central differences for x = 0(.01)3. The 
auxiliary function So(y) defined by the relation —oo(x) = xInx — x + So(y), 
= 1/x is given to 6D for y = 0(.01).35 with first differences. 

Comparison with the tables of the same functions published by NBSCL! 
showed a transposition in the value for x = .30 in the table of u(x) at the 
top of page 1. There is also a gradually increasing discrepancy of one to 

four units in the values of oo(x) from x = 2.40 to x = 3.00. 


MILTON ABRAMOWITZ 
NBSCL 


1NBSCL, Tables of Coulomb Wave Functions, v. 1. AMS no. 17, Washington, 1952 
[MTAC, v. 7, p. 101-102]. 


1162[L].—L. J. SLtaTer, “On the evaluation of the confluent hypergeo- 
metric function,’’ Cambridge Phil. Soc., Proc., v. 49, 1953, p. 612-622. 


This paper contains an 8S table of 


_ 
= + n)n! 


for a = —1(.1)1, 6 = .1(.1)1, x = 1(1)10, and the asymptotic expansion 
of this function for large x with converging factors which improve the 
accuracy quite considerably. 

The tables were computed from the power-series expansions, about 
thirty terms of the series being required. The computations were carried 
out on the EDSAC at the Cambridge Mathematical Laboratory. 


A. E. 


1163[L ].—UnIversity oF CALIFORNIA, Department of Meteorology, Tables 
relating to Rayleigh scattering of light in the atmosphere. Computed for 
the project Investigation of polarization of sky light by direction of 
ZDENEK SEKERA, Project Director. Computations performed under 
supervision of GERTRUDE BLANCH, by the U. S. Department of Com- 
merce, NBS, Los Angeles, INA. November, 1952, xxvi + 85 p., 1 p. 
Errata. 21 X 27 cm. 


This is a collection of nine tables useful in the study of the polarization 
of the light in the earth’s atmosphere. The numerical work presented here 
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is based on the theory developed by S. CHANDRASEKHAR,' who in 1947 
found the exact solution to this classical problem. A helpful feature of the 
publication is the careful definition of all tabulated quantities as well as a 
description of their physical significance. The fundamental functions on 
which the solution of the whole problem depends, are the four pairs of 
functions X®(u,7) and Y™(u, 7) (k = 1, 2, 3, 4) defined by the simul- 
taneous integral equations 

XQ) = 144 — a 
(1) 


where 
w(t) = +2 — 24), = (1+ 2° +244), 


v(t) = and y(t) = #(1 — 


Successive approximations X, (uz, 7) and Y,(u, 7) to the solutions of (1) 
were obtained by an iteration procedure, starting with the values 


Xo™(u, 7) = 1, Yo™(u, 7) = er. 


All calculations were made with IBM machines. Depending on the value 
of 7, six or seven iterations were carried out. Table I contains the values 
of X,(u, 7) and Y,(u, 7) to 5D, 

for r = 0.15 and 0.25, m = 6, w = 0(.02).32(.04).96(.02)1 

and for r = 1, m = 7, w = 0(.02).2(.04).96(.02)1. 
The intervals used in the numerical integrations were such that the entries 
are guaranteed only to within 0.0001 for » = 0.02, and to within 0.00003 
for »> 0.04. In order to exhibit the behavior of successive approximations 
to X and Y, in Table IX representative values are listed to 7D. Since for 


k = 3 the solutions of (1) are not unique, the functions tabulated (called 
“standard solutions’’) are those defined by the relations 


X® =X + ql X®% + YO]; Yo = y@ — quLx® + 
where 
x; and y;™ being the “moments” of the X and Y functions: 


1 
xy = f ve (t, dt, 


1 
= f v(t) (t, r)t? dt. 


Using the various formulas given by the theory, the authors of this 
publication finally illustrate, in Table VIII, the predicted behavior of the 
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sky radiation by considering the intensity and polarization on the principal 
meridian (containing the sun) for various elevations of the sun. The correc- 
tions to the intensity and polarization predicted by the theory for a ground 
surface reflecting according to Lambert’s law with albedos A = 0.10, 0.25, 
0.50, and 0.80, have also been included in Table VIII. 

The tables presented here not only provide the numerical results pre- 
dicted by an exact theory of atmospheric scattering according to Rayleigh’s 
laws, but also provide a basis for a quantitative evaluation of the non- 
molecular component of the sky radiation. 


G. Mtncu 
California Institute of Technology 
Pasadena, California 


1S. CHANDRASEKHAR, Radiative Transfer. Oxford, 1950. 
1164[L].—K.aus ZwEILinG, Grundlagen einer Theorie der biharmonischen 


Polynome. Verlag Technik, Berlin, 1952, viii + 130 p., 5 plates, 1 insert, 
18 X 24 cm. 


The polynomials tabulated in this volume are 
Poo = Pio =X, = xy 
= + n= 1(1)6, w= x+ ty, D= x — ty 


x(w" — w*) 


= 4ni 


nm = 1(1)5 


= + [2nx(w + — iy(w™ — nm = 1(1)5 


7 4n(2n + 2) 


Ponti = n = 1(1)5 


Pur = Pui, Peo = x Poi = x*y 


all for x, y = 0(.1)1. 
The zero lines of these polynomials are also tabulated, graphs of the 
P,,; are given, as are explicit representations in Cartesian and polar coordi- 


nates, various formulas and properties, and the expression of f PnP in dx 


in terms of products of these polynomials, for m + n,i+k< 6. 

The polynomials tabulated here are useful in problems concerning the 
bending of elastic plates. Similar quantities have also been tabulated by 
THORNE. 


A. E. 


1C. J. THornE, “A table of harmonic and biharmonic polynomials and their deriva- 


tives,” Utah Engineering Experiment Station, Bulletin, no. 39, 1949, supplement. Salt 
Lake City. 
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MATHEMATICAL TABLES—ERRATA 


In this issue references have been made to errata in RMT 1153, 1161, 
and in Nore 158. 


233.—Lr CENTRE NATIONAL D’EtuDEs DES TELECOMMUNICATIONS, Tables 
des fonctions de Legendre associées. Paris, 1952. [MTAC, v. 7, p. 178]. 


The authors report that on p. 17 the entry for m = 3.1, @ = 57° should 
read —0.4259087. 
A. E. 


234.—L. S. KHRENOV, Semiznachnyte Tabliisy Trigonometricheskikh Funkisit. 
(Seven-figure tables of trigonometric functions.) 1951 [MTAC, v. 7, 
p. 238-239]. 


In Table II, p. 131-401, the following errata were discovered by inspec- 
tion and checking of differences. 


Page Angle Function For Read 
153 3 34 50 csc 01284 01234 
176 7 26 30 tan 6160 6169 
196 10 47 0 sec 78757 79757 
202 11 41 40 ctg 1181 1179 
205 12 15 10 tan 1711 1721 
20 tan 2218 2228 
30 tan 2726 2736 
40 tan 3234 3244 
209 59 30 csc 835 935 
210 13 2 20 csc 828 928 
219 14 35 10 ctg 3878 2878 
220 44 10 csc 726 724 
250 19 45 30 sin 466 456 
251 51 20 sin 466 456 
277 24 15 50 ctg 8461 8481 
278 23 10 sec 9574 9554 
287 25 52 20 ctg 256 254 
300 28 7 50 sin 527 427 
313 30 11 40 ctg 86553 85553 
320 31 29 20 ctg 1778 1776 
336 34 3 40 csc 54683 54693 
343 35 12 40 ctg 1.4160069 0069 
50 ctg 8611 1.4168611 
18 20 cos 279 281 
347 58 10 csc 1136 1138 
376 40 43 20 ctg 6966 6956 
382 41 40 30 cos 9224 9284 


Twenty values in which the errors amounted to one unit in the last 
place have not been included in this list. 
S. A. JOFFE 
515 W. 110 Street 
New York 25, N. Y. 
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235.—NBS Applied Mathematics Series, No. 6, Tables of the Binomial 
Probability Distribution. Issued January 1950; reprinted October 1952. 


The following errors occur in one or both of the printings: errors marked 


with an asterisk occur in both editions; those without an asterisk were 
corrected in the reprinted edition. 


Entry 

Page n r p for read 
= 7 8 3 .06 .0086873 .0088773 
18 14 2 18 .2624913 .2724913 
116 37 0 .09 .0395163 -0305163 
134 40 3 .07 .2211640 -2311640 
134 40 4 .07 -1709448 -1609448 
140 41 5 .06 .0668162 .0628162 
140 41 6 .06 .0200573 .0240573 
192 49 17 a .0289183 -1189183 
192 49 18 .2038364 -1138364 
200 6 3 47 .6984534 .5984534 
*212 13 3,4, 5, 6; interchange values for p= .14and p=.15 
244 24 22 45 -0000001 -0000021 
*319 39 19 49 .5_ 68882 .5768882 
327 40 4 .07 .3163132 .3063132 
384 49 18 35 .5424066 .4524066 


It should be noted that several corrigenda may have resulted from a 
single error in computation or transcription. For example, the transposition 
of two digits on page 384 generated two further errors on page 192. 

We are indebted to many users of the tables for reporting these errors. 


CHURCHILL EISENHART 
National Bureau of Standards 
Washington 25, D. C. 
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177[A ].—J. W. WRENCH, JR., A New Approximation to 180 x-'. One manu- 
script page on deposit in the UMT FILE. 
This 2035D approximation is the by-product of the calculation of x. 
[See Note 159 in this issue. ] 


J. W. WRENCH, JR. 
4711 Davenport St., N.W. 
Washington 16, D. C. 


178[ F ].—A. FERRIER, Factor table for 3n* — 1. One photostat page 43 X 63 
cm. Deposited in the UMT FILe. 


This table gives 316 complete factorizations of 3n* — 1 for m < 1000. 
There is also a table giving the values of m modulo p for which 3n‘ — 1 is 
divisible by p = 12k + 1 < 3000. 

A. FERRIER 


College de Cusset 
Allier, France 
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179[L ].—T. H. Crow ey, Tables of Integrals of certain Bessel Functions. 
Available at the Antenna Laboratory, Ohio State University, Columbus, 
Ohio. 

These tables give values of the integrals 


wt cos 2 de 


for 
u = 0(.02)10, = 0(.1)10, Au< 15. 


Although the calculations were designed to give 4D accuracy, spot 
checking indicates an accuracy of 5D. 
T. H. CROWLEY 
Ohio State Univ. 
Columbus, Ohio 


180[L].—E. W. PrKe, Table of Parameters for the Summation Analogue 
of Laguerre Polynomials. Two typewritten pages on deposit in the UMT 
FILE. 


This table is of use in the design of filters for pulsed information. 


E. W. PIKE 
Raytheon Mfg. Co. 
Waltham, Mass. 
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TECHNICAL DEVELOPMENTS 
INPUT AND ORGANIZATION OF SUB-ROUTINES FOR FERUT 


1. Introduction.—Descriptions of methods for handling sub-routines on 
other machines have been written for the Manchester Electronic Computer,! 
EDSAC,? SEAC? and ILLIAC‘. The purpose of this article is to discuss the 
input and storage of routines and numerical data and the organization of 
routines during the solution of a problem on Ferut. 

Ferut is the serial, one address, electronic digital computer now in opera- 
tion in the Computation Centre, University of Toronto, Canada. It was 
built by Ferranti Ltd., Manchester, England, and is similar to the computer 
at the University of Manchester, England. Ordinary 5-hole telegraphic 
punched tape is used for input with a photoelectric reader. The output 
mechanism is a teleprinter and punch. 

Information in the machine is kept in blocks or pages. In each page 
there are 64 short lines of 20 binary digits each. For convenience these 20 
bits are arranged in four groups of five bits. The five-bit group constitutes 
a digit in the scale of 32 and is represented by one teleprint character. An 
electronic instruction consists of one short line, two teleprint characters for 
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the address, and two for the order. Numbers consist of two short lines or 
40 bits. 

The electronic storage comprises eight pages labelled So, Si, ---, Sz. In 
addition, there is a magnetic storage of 256 tracks on a drum, each track 
holding two pages of information. The pages may be transferred singly or 
in pairs to or from a magnetic track. This is accomplished by a magnetic 
transfer instruction. 

2. Routine Changing.—During the operation of a program, space in the 
electronic store is at a premium and thus the instructions necessary to carry 
out a problem are divided, as is customary, into a main, or master routine 
and a number of sub-routines. Only the set of instructions immediately in 
use is kept in the electronic store; the other routines are in the magnetic 
store. When it is desired to change from one routine to the next, it is neces- 
sary to have a group of orders which will (1) preserve, in a list called the link 
list, the location in magnetics of the routine being left, (2) record the line 
at which it is being left, (3) locate the next routine and call it down from 
magnetics and finally (4) enter the new routine at the appropriate line. 
This group of instructions is common to all routine changes and thus is kept 
separate from any one routine in a Routine Changing Sequence (R.C.S.) 
which is permanently down in the electronic store. (This occupies page S. 
along with powers of 2 for shifting and some other constants which are 
frequently required, and the page is known as PERM.) 

There must be, however, in the routine being left a few instructions 
supplying the required information about the specific routine to be entered 
and a control transfer to R.C.S. This linking sequence of orders which 
appears in the routine being left must be kept to a minimum to conserve 
space, and in the case of the Ferut system it occupies four short lines. One 
of these lines is, of course, the magnetic transfer which calls down the next 
routine from the magnetic store to the electronic store. 

This transfer is recorded in the link list by R.C.S. so that it can be used 
at a later date. The list of links grows as the control moves down a chain 
of sub-routines changing from the master routine (zero-th level) to the first 
level sub-routine, to the second level sub-routine, etc. The list is consulted 
by R.C.S. as the control returns back up the chain of sub-routines to the 
master routine. 

3. Directory.—The magnetic tracks are numbered, and in writing a 
magnetic transfer it is necessary to specify the track number. For several 
reasons it may be desirable to change tracks after a problem has been coded 
and the tapes have been punched. You must then seek out and alter all 
magnetic transfers which, of course, are embedded in the program. This is 
not only tedious but leads to errors. To overcome this difficulty all magnetic 
transfers are kept separate from the routines proper in a list called the 
directory. In the line of the routine where the transfer should appear is 
placed the address in the directory of the required magnetic transfer. This 
line is called a false line. During input the false lines are replaced by the 
appropriate magnetic transfers from the directory. 

If the magnetic transfer required to call down a certain routine is stored 
in line ab of the directory, ab is said to be the Directory Number of the 
routine. This directory number is what is written in a false line. When a 
routine has been tested and is considered useful, it may be assigned a 
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Library Number which will be its permanent directory number. Sixteen lines 
at the beginning of the directory are kept open for routines of a particular 
problem. These directory numbers are temporary. 

4. Input and Output.—To perform input or output operations it is nec- 
essary to combine the simple input or output orders of the machine with 
other orders in an input or output routine. The function of the Ferut Input 
Routine is to read information from tape, perform certain alterations on 
routines or numerical data and store the routines or data in assigned loca- 
tions in the machine. 

The input routine reads one character at a time from the tape testing 
it until it finds a Warning Character. There are eleven warning characters. 
Each of the warning chs. causes control to be taken to a certain section of 
the input routine where a prescribed group of instructions is obeyed. 

An attempt has been made to choose, for warning chs., teleprint char- 
acters which are suggestive of the operations which the warning chs. perform. 
For example—the warning ch. J followed by n/ causes the n chs. following 
the / to be printed out by the machine. (The ch. J suggests the word JOT.) 
After the Jn/ the machine reads in n chs. before it begins again to examine 
the tape for the next warning ch. Since / (zero) is not a warning ch., blank 
tape is examined rapidly and passed by. Thus between meaningful sequences 
on the tape (such as the Jn/ followed by n chs.) it does no harm to have 
spaces, and they are a great assistance when reading the tape by eye. 

The use of the electronic pages during input is best described by a 


diagram. 
So Si Ss Ss Ss Ss 
INPUT ROUTINE PERM LINK RECEIVES DIRECTORY 
LIST INPU 


NPUT 
INFORMATION 


After input only page S2 and one half of S; are used. Routines normally 
operate from So and S;. 
The list of warning characters is as follows: 


STOP—Stops the tape—control is returned to input to search for the next warning ch. 
when a prepulse is given manually. 
Jn/—Prints out the n chs. on the tape following Jn/—used for titling routines 
(J suggests JOT). 
E ab—Transfers control to line (ab + 1) (E suggests ENTER). 
: abcd—Causes abcd to be obeyed as a magnetic instruction (: suggests /: the “obey 
as a magnetic instruction’”’ order of Ferut). ; 
K ab n—Reads in from the tape n short lines starting the entries in line ab of the store. 
A maximum of 32 short lines can be read in at a time. (K suggests absolutely 
nothing but is used for this purpose by all the other input systems.) 


sao sl the routine whose directory number is ab from 


pages S, and §S; if it is a 2 page routine) oe the magnetic track assigned 
in the directory. Check of the magnetic transfer is performed, and a warning 
signal indicates failure. 


page S, (or 
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cats in as an the routine whose directory number is ab and 
enters it at line (cd + 1). (@ sounds like ad—and R is suggestive of Routine.) 
F ab—Replaces the false line ab of a routine by the contents of a line of the directory. 
(Explanation: In order to incorporate a magnetic transfer within a routine in 
line ab the directory number of the appropriate magnetic transfer is placed in 
the routine itself in line ab. This is a false line. The warning ch. F replaces 
this by the appropriate true line, i.e., the magnetic transfer as it appears in 
the directory. For example: One line of the linking sequence is such a false 
line and any routine containing a linking sequence must be followed by the 

proper F sequence.) 
C/(2n)—Adds together the 32 long lines of page n and compares the sum modulo 2@ 
with the known Check Sum of this page. If the two differ, a signal occurs. The 
known check sum must, of course, be planted in a certain line of the store. 


A wide variety of routines have been devised for handling the input and 
output of numerical data depending on the position of decimal or binary 
point, method of conversion, number of significant figures, layout, etc. 
These are entered as sub-routines from the input routine by an R warning 
ch. or from the master routine of the problem by a linking sequence. When 
numerical data are put in, the pages of information can be given directory 
numbers and thus stored and brought down in a manner similar to that 
used in handling sub-routines. 

The following checks are a part of this system of input. Each Up or 
Down transfer is checked; thus the contents of magnetic and electronic 
storage are known to be identical. The transfer from tape to electronic 
storage is checked by having a known check sum of a routine or of numerical 
data read in and compared with the check sum of the contents of the elec- 
tronic page by means of the C warning ch. Alternately, if the check sum is 
not known, the tape may be read in twice, once with current on which 
writes from electronic to magnetic storage, the second time with this write 
current off. The Up transfer check then effectively compares what was read 
in the first time from tape with what was read in the second time. 

Most problem tapes consist of sequences of characters which (i) print 
the title, (ii) read in and store the directory, (iii) read in, alter and store 
routines and (iv) start the problem. Titling (see i) is useful for identifying 
tapes. Process (ii) consists of reading the magnetic transfers for the routines 
being used into the appropriate lines of the directory. The directory is 
stored by means of the : warning ch. All routines are read (see iii) by K- 
sequences into S, (or S, and S; in the case of two-page routines). They are 
checked if desired by the C ch., altered if necessary, i.e., false lines are 
replaced by F-sequences, and then stored by the appropriate Up order. 
The problem is started (see iv) by entering the master routine as an ad- 
routine by means of an @ sequence on the tape. 

Routines may be altered during input by bringing them down with a 
D-sequence, altering the desired lines by K-sequences and then sending 
them back up to magnetics by a U-sequence. 

Library tapes may be used without alteration in a problem tape. They 
consist of a titling sequence, K-sequences reading the routine into S, or S, 
and S; and then the Up order using the library number of the routine. 
Library tapes have a check sum taken by means of the Down order followed 
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by the C warning ch. After this false lines are replaced by F-sequences, and 
the routine is stored again by the U-sequence. 

5. Example.—Consider the organization of the routines for the two- 
dimensional Fourier synthesis 


o(y, 2) = Fucos 2x(ky + Is) 
k=0 l=0 


where the coefficients F,, are tabulated. The program for this problem will 
consist of four routines and the data. Three of the routines to be used are 
library routines, namely, Decimal Input, Decimal Output, and Cosine. 
The fourth is the Master routine. The latter must call in Cosine and Decimal 
Output as sub-routines and thus has a false line in each of the linking 
sequences calling in these two routines. These false lines have addresses 
ayb,; and abe. Transfers calling down the data require false lines in a;b;, 
aab,, etc. 

The coefficients F,; are taken in and converted from decimal to binary 
form by the Decimal Input routine before they are stored in the magnetics. 
This is accomplished during input by entering Decimal Input as a sub- 
routine from the tape by means of the R warning ch. Details for the routines 
and data are given in Table I. 

The magnetic transfer calling down the Master will be (t:) 52So, where 
(t:) represents the two teleprint character track address, So indicates the 
electronic storage from which the routine operates, and 62 denotes that it is 
a two-page routine. This magnetic transfer must be recorded in line /N of 
the directory. 


TABLE I 
Directory Magnetic Electronic 
Description Number Storage Storage 
Master /N ti So, Si 
Decimal Input SF te So, S1 
Decimal Output ON ts So, Si 
Cosine XN t, (left half) So 
Data 1 (wo pages} EN ts Ss, Ss 
Data 2 (two pages @N ts Se, Sz 
Data 3 (two pages) AN tr Sa, Ss 
etc. 
The tape is organized in the following way: 
Punching Description 
J 33 3 T is a letter which is not a warning ch. but 
indicates the beginning of the tape to the 
operator. 
JZ / FOURIER: Title (Z is the teleprint character for 17.) 
SYNTHESIS 
The K-warning ch. reads into the directory 
the magnetic transfer for bringing down: 
K / NE (t1) 62 So the Master 
K SFE (te) 82 So Decimal Input 
KONE (ts) &2 So Decimal Output 
KXNE (t,) 6: So Cosine 
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Punching Description 
KENE (ts) 83S, Data 1 
K@NE (te) 52 Ss Data 2 
KANE (tz) 62S. Data 3 
etc. (E is teleprint for 1.) 
:I1/RN I /RN is the magnetic transfer instruction 
which stores the directory on track. 
STOP 
Library tape for Decimal Input 
Library tape for Decimal Output these routines all have check sums 


Library tape for Cosine 
Master routine punched with the K-sequences which read it into S, and Ss. 


Faib; Fasbz Fasbs Faby, etc. Replaces false lines in Master by the mag- 
netic transfers from the directory. 
U/N Stores Master on track ty. 
STOP 
RSFcd Enters Decimal Input as a sub-routine at line 
Data 1 cd +1 whereupon Data 1 punched in the 
appropriate form are read in and converted. 
UEN Stores Data 1 on track ts. 
RSFcd Reads in and converts Data 2. 
Data 2 
U@N Stores Data 2. 
RSFcd Reads in and converts Data 3. 
Data 3 
UAN Stores Data 3. 
etc. 
STOP 
@/Ncd Enters the Master as an ad-routine at line 
ed + 1. 
$$S$$.... A non-warning ch. indicating the end of the 
tape. 


6. Summary.—The Ferut System has the following features: 

(1) No reference to fixed magnetic tracks is made except in a special 
list called the directory. 

(2) Minimum time and mental effort are consumed in changing from 
one routine to another during the operation of a problem. 

(3) The details of the input are easy to remember. 

(4) Alterations may be made to routines or pages of data with a mini- 
mum of effort. 

(5) Certain automatic checks are built into the system. 

(6) Library tapes may be used without alteration in a problem tape. 

This system of input and routine organization was adopted after study- 
ing the various systems in use at the University of Manchester.' The details 
were devised by the author in cooperation with C. Stracuey, National 
Research Development Corporation, London, England, and H. GELLMAN, 
Chalk River, Canada, with the advice and criticism of C. C. GotLieB and 
the staff of the Computation Centre, University of Toronto. 


J. N. P. Hume 


University of Toronto 
Toronto, Canada 
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1 The Programmers Handbook for the Manchester Electronic Computer (Mark II). mimeo- 
graphed, Univ. of Manchester, England, Aug. 1952, 9 chapters. 
2M. V. Wirkgs, D. i WHEELER & S. . The Preparation of Programs for an Elec- 
tronic Digital Computer. ddison-Wesley Press, Inc., Cambridge, Mass., 1951. 
D. R. HaRTREE, Numerical Analysis. Chapter XI, Oxford University Press, 1952. 
3 Anon., “The incorporation of subroutines into a complete problem on the NBS Eastern 
Automatic Computer,” MTAC, v. 4, 1950, p. 164-168. 
*D. J. WHEELER, Program Organization for the University of Illinois Digital Computer. 
Internal Report No. 33, University of Illinois, May 16, 1952. 


DISCUSSIONS 
Programs for Computing the Hypergeometric Series 
The hypergeometric series 


abs | a(a+1)b(6+1) 
c 1! c(c + 1) 2! 


appears in the solution of various problems in applied mathematics and it is 
therefore desirable to have efficient methods for evaluating this series. To 
meet these needs a program of investigating available methods for com- 
puting the hypergeometric series is being conducted. The present note 
discusses two computations on SEAC. 

A code for real values of parameters and argument has been checked 
out with a, b, c equal to various combinations of 1,2,3 and —.9 SxS .9 
where many of the answers are known functions. The results were accurate 
to nine and more often ten significant figures. For each set of parameters 
it took approximately twenty-five minutes to tabulate the series for 
x = —.9(.1).9. 

A code for complex values of parameters and argument has been checked 
with one of the real cases above giving the same accuracy, and with 
a@=1+%2.2,b = -2,c =2+%, ¢=.5 —4.8 where the hypergeometric 
series is a polynomial. The agreement in this case was exact. 

A previous hand computation of F(a, b; c; z) for 


@ = —11.753 + 712.204, b = 12.753 — 412.204, c=2, 2=.5 
gave the result 


F(a, b3¢;2) =1+— 


F(a, b;c;z) = 832,109 — 1827,535 


where the result was accurate to four figures. Using these same parameters 
and argument and carrying all computations to at least five decimal places 
and computing until the truncation error was less than 10? (54 terms) 
SEAC gave 

F(a, 6; c; z) = 832,103.8 — 1827,536.8. 


Repeating the computations to at least seven decimal places and computing 
first until the truncation error was less than 10-* (57 terms) and then until 
the truncation error was less than 10-* (65 terms) SEAC gave the identical 
results 

F(a, c; 2) = 832,108.0985 — 71827,540.7310. 


As a check on this result F(a — 1,6;c;2) and F(a,b + 1;c + 1;2) 
were computed and the following recursion formula yielded the results 
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F(a, b;c;2z) = F(a — 1,6;c;2) + bzF(a, b + 1;¢ + 1;2)/c 
F(a, 6; c; 2) = 827,563.5390 — 146,025.5154 

+ 62(123,187.2697 — 4127,239.3492)/c 
F(a, 6; c; z) = 832,108.0972 — 1827,540.7303 


which agrees with the previous computation to the second decimal place. 
Each of the above computations took four minutes time on SEAC, including 
reading in and out. 

This code uses double precision operations and begins by reading in 
a, b, c, z and the size of the term where the series is truncated, where all are 
twenty-two digit decimal numbers. After computing and printing the result 
the SEAC calls for a new z and will compute using the previous parameters, 
or may be made to read in new parameters by sending the Control to address 
zero. In general the result will be accurate to 10"~§ significant figures where 
F(a, b; c; z) and its largest term are both less than 108. 

This last code has been modified to compute the confluent hypergeo- 
metric function for complex parameters and argument. Using the modified 
code w = 1.35218 220994 has been found as a zero of the function’ 


1-—w 
:1;20), 


NBSCL 


1H. A. LAUWERIER, “The use of confluent hypergeometric functions in mathematical 
physics and the solution of the eigenvalue problem,” Applied Scientific Research, 1951, A, 
z, p. 184-204. 


W. F. 
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1091. Anon., “Analog-to-Digital Converter simplifies data reduction,” 
Product Engineering, v. 24, May 1953, p. 234-235. 


Principles of the SADIC analog-to-digital converter built by Consoli- 
dated Engineering Corp. are briefly described. 
R. D. ELBouRN 
NBSCML 


1092. ANon., ‘““Computer assembly line,” Electronics, v. 26, July 1953, p. 
200. 


This is a short expository paragraph about the IBM electronic data- 
processer Model 701. 


1093. ANon., “Converters put data into useful form,” Aviation Week, v. 58, 
May 18, 1953, p. 68. 


Analog-to-digital converters having three decimal digit precision and a 
cycle slightly less than one second have been announced by two companies, 
the SADIC by Consolidated Engineering Corp. and the Teleducer by Tele- 
computing Corp. 

R. D. ELBouRN 
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1094. ANON., “Mathematicians needed to help ‘brain,’ Machine Design, 
May 1953, v. 25, p. 236. 
Anon., “Armour Research Foundation to get electronic super-brain,”’ 
Midwest Engineer, 5, Feb. 1953, p. 28. 


Both articles deal with the IBM electronic computer 701. The first gives 
a brief description of this card-programmed calculator with its magnificent 
equipment of electrostatic memory, magnetic drum, and magnetic tapes. 
The second announces the prospective rental of this model by the ARF of 
the Illinois Institute of Technology. 

The pleasure we find in reading of any advancement in the field of elec- 
tronic computation is considerably dampened by the persistent occurrence 
in print of such epithets as “super-brain,”’ “thinking machine,” and the 
like—as applied to computers. Nor is it a pleasure for the hard-working 
mathematician to be told that ‘fantastically hard problems are easily solved 
by the electronic computers.” 

It should be realized that man-made machines cannot solve problems, 
any more than telephones can carry on conversations, or telescopes can 
make celestial observations. None but the human mind is capable of initi- 
ating and solving mathematical problems, as well as of utilizing various 
types of equipment for the computation of the results. 

Our ruffled feathers were smoothed somewhat by the timely warning 
voiced in the first article that “‘there may not be enough top quality mathe- 
maticians to prepare problems for . . . these machines. . . .”” Prophetic 
and true words! 


I. RHODES 
NBS 


1095. ANon., “NBS designs ‘building block’ computer,’’ Product Engi- 
neering, v. 24, March 1953, p. 205. 


A new computer, under development at the National Bureau of Stand- 
ards for the Defense Department, is composed largely of approximately 800 
standardized circuit packages. These are expected to take care of 90 percent 
of the total circuitry, exclusive of memory units. Two general types of 
standardized packages have been developed ; the first contains an amplifier 
tube, pulse transformer and germanium diodes; and the second, delay lines. 
These basic types are made in several versions as dictated by circuit re- 
quirements. 

The assemblies, now being manufactured commercially, appear to be 
about two inches wide, five inches high and 10 inches deep. Connection is 
through a 60-contact pin assembly at the rear. To facilitate checking and 
fault-localization, a test jack is located on the front of each unit. 


M. L. GREENOUGH 
NBS Electronic Instrumentation 


1096. AN WANG, “High-speed number generator uses magnetic memory 
matrices,” Electronics, v. 26, May 1953, p. 200-204. 


Numbers are displayed in an 8X8 dot pattern on a cathode ray tube 
by using conventional raster-scan deflection and pulsing the beam on 
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according to character shapes stored in 8X8 magnetic core matrices. The 
device has been tested at 8,000 characters per second, but the builder 
believes 100,000 per second is achievable. 


R. D. ELBourn 
NBSCML 


1097. W. Cocuran & A. S. Douctas, “A new application of EDSAC to 
crystal structure analysis,’ Nature, v. 171, June 20, 1953, p. 1112-1113. 


Determination of the structure of centro-symmetric crystals by evalu- 
ating the electron density for each of the 2” possible sign combinations on 
N structure factors can be made after the correct signs of about twenty of 
the largest structure factors have been found (for crystals of moderate 
complexity). 

The EDSAC surveys this large number of alternative answers and selects 
from them a few which fulfill certain criteria. In the case of salicylic acid, 
it completed in less than three hours the calculation for 2° = 65,536 inde- 
pendent possibilities. The number of sign combinations corresponding to 
each such possibility is listed in an accompanying table in this article. This 
method is considered valuable when used in conjunction with PEPinsky’s 
X.R.A.C. for rapid testing of the most plausible sign combinations. 

Mention is made of another method using EDSAC which requires only 
a few minutes of machine time for determining rapidly the structures of 
crystals of moderate complexity. 


E. MARDEN 
NBSCL 


1098. P. Kiass, ‘“‘Giant brains could aid air defense,” Aviation Week, v. 58, 
May 11, 1953, p. 67-68+. 


This article is a description of two new computing machines, Reming- 
ton Rand’s 1103 and International Business Machines Corporation’s 701 ; 
it stresses their possible uses in the aviation field. The two uses mentioned 
were both ‘“‘real time” control uses, i.e., air defense operations and airport 
traffic control. No mention was made of their main use so far, that of solving 
complex problems and simulating flight conditions. ‘Real time’’ control can 
consist of either giving fast answers pictorially and numerically to human 
beings or of actually controlling aircraft and guided missiles by remote 
control. Aside from the vast amount of pre-stored information, information 
as to changing conditions can be fed into the computer directly from radar 
and other equipment, thus enabling the machine to make its lightning 
computations. 

The bulk of the article contrasts the specifications for the two machines, 
giving a good summary of their technical points. The machines are roughly 
comparable but quite different in some details. 

C. J. Swirt 
NBSCL 
NEws 


Symposium on Monte Carlo Methods.—A Symposium on this subject is being spon- 
sored by’ the AeronauticalJResearch Laboratory, WADC. It will be held at the University 
of Florida at Gainesville, on March 16 and 17. Those interested in Monte Carlo methods 
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are invited to attend. Further information may be obtained from H. A. MEvER, University 


of Florida, Gainesville. 


Association for Computing Machinery.—A general meeting of the Association was held 
at the Massachusetts Institute of Technology on September 9, 10, and 11, 1953. Many 
new ideas and developments in the applications of computing machinery were reported in 
parallel sessions. In addition the Whirlwind I Computer located at MIT was demonstrated 
for those attending. The program for the meeting was as follows: 


September 9, 1953, Morning Session 


Punched card techniques 

Solution of simultaneous linear equations by 
the Crout method on the IBM card pro- 
grammed calculator 

Remington Rand 409 Electronic Computer 
as applied to scientific and commercial 
problems 

Matrix multiplication on standard punched- 
card machines 

A calculation system for the IBM CPC 
Model 1A 


The EIP—An external and internal pro- 
gram setup for IBM’s Model II CPC 

A three factor floating decimal panel for the 
IBM CPC 

A high-speed multi-purpose board for the 
IBM card-programmed electronic calcu- 
lator 

Planned “work flow” in an engineering com- 
puting installation 


F. M. Verzun, MIT, Chairman 
K. J. Berc, United Aircraft Corp., East 
Hartford, Conn. 


R. T. Bruce, Remington Rand, Inc., 
South Norwalk, Conn. 


James Morison, Douglas Aircraft Co., 
Inc., Santa Monica, Calif. 

J. C. Rea, Engineering Calculations 
Group, Allison Division, General Mo- 
tors Corp., Indianapolis, Indiana 

J. C. Smaw, The Rand Corp., Santa 
Monica, Calif. 

D. S. Byron, Bell Aircraft Corp., Buffalo, 
N. Y. 

J. A. PostLEy, Hughes Research and De- 
velopment Laboratories, Culver City, 
Calif. 

Rex Rice, Jr., Northrup Aircraft, Inc., 
Hawthorne, Calif. 


September 9, 1953, Morning Session 


Digital computer techniques 
Single vs. triple address computing machines 


A method of radix conversion 
Interim and comprehensive systems of com- 
putation on MIDAC 


Matrix multiplication on the ERA 1103 


Digital computer mathematics 


A new method of determination of all roots 
of an algebraic equation 

A comparison of machine methods for eval- 
uating certain mathematical functions 


The numerical solution of semilinear para- 
bolic equations by difference methods 


C. W. Apams, MIT, Chairman 

C. C. Extcot, Naval Ordnance Labora- 
tory, White Oak, Md. 

WALTER Sopen, Analysis Division, U. S. 
Naval Air Missile Test Center, Point 
Mugu, Calif. 

J. H. Brown and J. W. Carr, Univ. of 
Mich., Willow Run Research Center, 
Ypsilanti, Mich. 

A. E. Roserts, Engineering Research 
Associates, Division of Remington 
Rand, Inc., Arlington, Virginia 

J. W. Remington Rand, Inc., 
Chairman 

J. W. Carr, Univ. of Mich., Willow Run 
Research Center, Ypsilanti, Mich. 

W. B. Fritz, Ballistic Research Labora- 
tories, Aberdeen Proving Ground, Aber- 
deen, Md. 

A. J. Perwis, Purdue Univ., Statistical 
Laboratory, Lafayette, Indiana 
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The continued-fraction algorithm for com- G. W. KiNG, International Telemeter 


puting machines Corp., Los Angeles, Calif., and L. J. 
Epson, Arthur D. Little, Inc., Cam- 
bridge, Mass. 

September 9, 1953, Afternoon Session 
Numerical analysis F. B. HitpEBRAND, MIT, Chairman 
A survey of methods of solving systems of J. H. Curtiss, New York Univ., New 
linear algebraic equations York City 
A survey of methods for dealing with eigen- | H. H. Gotpstine, The Institute for Ad- 
value problems vanced Study, Princeton, N. J. 
Problems on partial differential equations F. J. Murray, Department of Mathe- 
matics, Columbia Univ., New York 
City 
September 9, 1953, Afternoon Session 
Recent systems developments C. C. Hurp, IBM Corp., Chairman 
MIDAC—The Michigan Digital Automatic J. E. DeTurx, R. Hock, J. Kaurman, & 
Computer H. BETHEL, Univ. of Michigan, Ann 


Arbor, Mich. 
The IBM Magnetic Drum Calculator Type _F. E. Hamirton and E. C. Kusie, IBM 
650 Corp., Endicott, N. Y. 
Storage class control in the ERA 1103 W. G. WELCHMAN, Engineering Research 
Associates, Div. of Remington Rand, 
Inc., Arlington, Va. 
The problems of preparing acceptance tests P. Brock & S. M. Rock, Consolidated 


for digital computers Engineering Corp., Pasadena, Calif. 
Automatic strain-gage and thermocouple re- R. PERLEy, United Aircraft Corp., East 
cording on punched cards Hartford, Conn. 
The analyzing reader D. H. SHEPARD, Intelligent Machines Re- 
search Corp., Arlington, Va. 
Computation for numerical control J. H. Runyon, Servomechanisms Lab- 


oratory, MIT, Cambridge, Mass. 
Large-scale computer output through the L. W. SprinKLE, Mathematical Compu- 
electronically controlled typewriter tation Branch, Air Force, Alexandria, 
Va. 


September 10, 1953, Morning Session 


Digital computer programming A. S. HousEHOLDER, Oak Ridge National 
Laboratory, Chairman 
An interpretive routine for the manipu- D. E. Mutter, Digital Computer Lab- 


lation of expansions in Boolean algebra oratory, Univ. of Illinois, Urbana, IIl. 

The IBM 701 speedcoding system J. W. Backus, New York Scientific Com- 
puting Service, IBM Corp., N. Y. 

The editing generator A. M. Koss & J. H. Warre, Remington 


Rand, Inc., Philadelphia, Pa. 
Analytical differentiation by a digital com- H. G. KAHRIMANIAN, Remington Rand, 


puter Inc., Philadelphia, Pa. 
Digital computer programming M. V. Witkes, Cambridge, England, 
Chairman 
Résumé of automatic coding techniques J. M. Bennett, Ferranti Ltd., Manches- 
being developed for digital computers ter, England 


N. Rocuester, IBM Corp., Poughkeep- 
sie, N. Y. 


| 


September 10, 1953, 


Punched card mathematics 

A method for solving boundary value prob- 
lems of mathematical physics on punch 
card machines 

A method of determining plate bending by 
use of a punched card machine 

Computing supersonic flow around axial 
symmetric bodies using method of char- 
acteristics 

Control panels for a card programmed 604 
computing system 

Logical algorithms 

Derivation of a prediction function for com- 
plex system performance 

Extension of the Veitch chart method in 
computer design 

Automatic series-parallel circuit simplifica- 
tion 

Automatic analytic geometry in machine 
design 


September 11, 1953, 


Numerical solution of partial differential equa- 
tions 
Finite-difference approximations to the fun- 
damental frequency of a vibrating mem- 
brane 
Boundary conditions in random walks 


Numerical treatment of a fourth order para- 
bolic partial differential equation 
Numerical solutions of some nonlinear heat- 
transfer problems 
Operation of a computation center 
Operation of a technical computing facility 


United Aircraft Corporation—Trinity Col- 
lege, Computation Course 


The educational aspects of operating a com- 
puting service 


AUTOMATIC COMPUTING MACHINERY 


G. M. Hopper, Remington Rand, Inc., 
Philadelphia, Pa. 

C. W. Apams, Digital Computer Labora- 
tory, MIT, Cambridge, Mass. 

S. Git, University Mathematical Lab- 
oratory, Cambridge, England 

D. J. WHEELER, Univ. of Illinois, Urbana, 
Til. 

J. W. Carr, III, Univ. of Michigan, 
Ypsilanti, Mich. 


Morning Session 


C. C. Lin, MIT, Chairman 

S. BercmaNn, Applied Mathematics and 
Statistics Laboratory, Stanford Univ., 
Palo Alto, Calif. 

A. D. WasEL, Mathematics Dept., Univ. 
of Santa Clara, Santa Clara, Calif. 

M. Rosrnson, Bell Aircraft Corp., Buf- 
falo, N. Y. 


N. A. Patton, Lewis Flight Propulsion 
Laboratory, Cleveland, Ohio 

S. H. MIT, Chairman 

E. D. FuLLENwiper, Missile Evaluation 
Laboratory, NBS, Corona, Calif. 

R. W. Brooks, Raytheon Manufacturing 
Co., Waltham, Mass. 

R. J. Netson, IBM Corp., Endicott, 
N. Y. 

L. P. Crosman, Laboratory of Advanced 
Research, Remington Rand, Inc., Nor- 
walk, Conn. 


Morning Session 


R. F. Crrpprncer, Raytheon Manufac- 
turing Co., Waltham, Mass., Chairman 

G. E. Forsyte and B. F. Hanpy, NBS, 
Los Angeles, Calif. 


W. F. Baver, Univ. of Michigan, Ypsi- 
lanti, Mich. 

S. H. Cranpati, Dept. of Mechanical 
Engineering, MIT, Cambridge, Mass. 

F. Att, NBS, Washington, D. C. 


W. J. Ecxert, IBM, Chairman 

W. D. BELL, Telecomputing Corp., Bur- 
bank, Calif. 

W. RamsHaw, United Aircraft Corp., 
Machine Computation Laboratory, 

East Hartford, Conn. 

G. T. Hunter and D. R. Mason, IBM 

Corp., N. Y. 
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Installation of an IBM 701 


H. R. J. Groscu, General Electric Co., 
Evendale, Ohio 


September 11, 1953, Morning Session 


Analog computation 


Automatic computation by analogue com- 
puters of certain errors in flight simulation 
problems 

Solution of certain partial differential equa- 
tions by analogue computers 

Analog computer solution of partial differ- 
ential equations 


Analog computer approach to heat transfer 
problem 

A method of approximating the real roots 
of a polynomial using an analogue com- 
puter 


Analog computation 
Equipment reliability as applied to analogue 
computers 


Generation of a function of two variables 


Generator for functions of two variables 


Survey of analogue multiplication 


An operational digital divider 


G. A. Philbrick Research 
Associates, Chairman 
G. L. Lanpsman, NBS, Corona, Calif. 


O. L. Bow1e, Watertown Arsenal, Water- 
town, Mass. 

H. B. Betck, Rensselaer Polytechnic In- 
stitute, Computer Laboratory, Troy, 
N. Y. 

C. F. Kayan and V. Pascuxts, Columbia 
Univ., N. Y. 

R. E. Carrot, Bell Aircraft Corp., Buf- 
falo, N. Y. 


W. W. Serrert, MIT, Chairman 

H. Jacoss, Jr., MIT, Dynamic Analysis 
and Control Laboratory, Cambridge, 
Mass. 

R. P. JERRARD and G. T. Jacos1, General 
Electric Co., Schenectady, N. Y. 

G. R. Wett1, MIT, Dynamic Analysis 
and Control Laboratory, Cambridge, 
Mass. 

C. M. Epwarps, Bendix Research Lab- 
oratory, Detroit, Mich. 

M. A. Meyer, B. M. Gorpon, and R. N. 
Nicota, Laboratory of Electronics, 
Inc., Boston, Mass. 


September 11, 1953, Afternoon Session 


Operation of a computation center 
The University Computing Laboratory 


Operation of the Office of Statistical Services 
Operating a computer efficiently 


Operating and administrative procedures for 
MIDAC 

Problems arising in the administration of a 
multi-machine digital computing service 

An integrated data processing facility 


C. V. L. Smita, ONR 

A. J. Pertis, Purdue Univ., The Univer- 
sity Computing Laboratory, Lafayette, 
Ind. 

F. M. Verzun, MIT, Office of Statistical 
Services, Cambridge, Mass. 

C. C. Gorires, Computation Centre, 
Univ. of Toronto, Toronto, Canada 
D. NEEB, Univ. of Michigan, Ypsilanti, 

Mich. 

J. W. Fiscupacu, Ballistic Research Lab- 
oratories, Aberdeen, Md. 

E. M. McCormick and H. H. Rosen, 
Missile Evaluation Laboratory, Na- 
tional Bureau of Standards, Corona, 
Calif. 
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September 11, 1953, Afternoon Session 


Business data handling R. A. Manetnt, John Hancock Life In- 
surance Co. 
Premium billing performed by large-scale R.T. WisEMAN, Sun Life Assurance Co. 
computers of Canada, Montreal, Canada 
Periodic billing and accounting E. F. Coo.ey, Prudential Insurance Co. 


of America, Newark, N. J. 
The development and application of elec- E. J. CunnrncHAM, Monsanto Chemical 


tronic equipment in Monsanto’s account- Co., St. Louis, Mo. 
ing department 
Experience with the Census UNIVAC D. H. HEtser & J. L. McPHERsoN, Bu- 


reau of the Census, Department of 
Commerce, Washington, D. C. 

Budget computation on IBM 701 Lt. R. E. Urman, Navy Aviation Supply 
Office, Philadelphia, Pa. 


OTHER AIDS TO COMPUTATION 
BIBLIOGRAPHY Z 


1099. ANnon., ‘Chart for the resolution of angles,” Product Engineering, 
v. 22, no. 6, 1951, p. 187. 


To be used for representing spatial figures on a plane. 


1100. Anon., ‘‘High speed response in electromechanical integrator,’’ Prod- 
uct Engineering, v. 24, no. 6, 1953, p. 240. 


A combination of variable speed drive and d.c. feed back integrator 


constitutes an integrator whose input is a voltage regarded as a function of 
the time and with mechanical counter output. 


1101. E. BromBerc & R. D. McCoy, “Calculating machines—new tools 
for the designer,” Product Engineering, v. 22, no. 3, 1951, p. 85-88. 


This article contains brief descriptions and photographs of REAC equip- 
ment and components including the servo multiplier and resolver and 


function generator. Also various applications to engineering problems are 
described. 


1102. F. W. Buss & W. L. Morris, “How analogical computing devices 
can serve process industries,” Chem. Engineering, v. 57, no. 7, 1950, 
p. 142-144. 


1103. W. H. Burrows, ‘Methods of calculating with graph papers,” Prod- 
uct Engineering, v. 22, no. 4, 1951, p. 140-145; also no. 6, p. 168-171. 


1104. F. P. Cozzone, “Organizing a computer center in the engineering 
department,”’ Product Engineering, v. 23, no. 1, 1952, p. 136-141. 


Various types of digital and analogue computers are described. A chart 
is given, showing the relative suitability of these types for aircraft design 
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problems. The organization, personnel problems, power requirements and 
space layout for an engineering computing center are considered. 
F. J. M. 


1105. A. S. Hatt & D. C. Tao, “Analysis of a symmetrical five-bar link- 
age,’ Product Engineering, v. 23, no. 1, 1952, p. 175-177; also “Design 
charts for a five-bar linkage,” p. 201, 203, 205. 


1106. B. H. List, R. C. McMaster, & R. L. MERRILL, ‘“‘Analogous systems 
in engineering design,”” Product Engineering, v. 24, no. 1, 1953, p. 184— 
195. 


Various network type analogues and a differential analyzer are dis- 
cussed. 


1107. C. P. Nacnop, “‘Nomograph for the volume of cones or pyramids,” 
Product Engineering, v. 22, no. 5, 1951, p. 209. 


1108. D. W. PEacemMAN & J. E. Vivian, ‘Bantam differential analyzer,” 
Chem. Engineering, v. 57, no. 8, 1950, p. 106-107. 


This article describes a small differential analyzer built by A. B. MACNEE 
and using his multiplying unit. 


1109. D. H. Pickens, “The electronic analog computer,” Product Engi- 
neering, v. 24, no. 5, 1953, p. 176-185. 


Principles and certain engineering applications of an electronic differ- 
ential analyzer based on a feed back d.c. amplifier are discussed. 


1110. W. W. Soroka, “Equivalent dynamical systems for complex vibra- 
tion problems,” Product Engineering, v. 23, no. 7, 1952, p. 130-133. 


A method of setting up a mechanical analogue with lumped masses and 
linear springs for structural vibration problems is described. 


1111. W. W. Soroka, “Resistance network analogue for solving vibration 
problems,” Product Engineering, v. 22, no. 4, 1951, p. 103-105. 


A manually adjusted resistance network is described. 


1112. E. C. VArnum, “Circular nomogram theory and construction tech- 
nique,’ Product Engineering, v. 22, no. 8, 1951, p. 152-156. 


An explanatory article with engineering examples. 


1113. E. C. VARNum, ‘“‘Nomogram for evaluating test data,” Product Engi- 
neering, v. 24, no. 2, 1953, p. 215. 


A circular nomogram for applying the “‘t test’’ to determine the signifi- 
cance of a difference of means. 
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NOTES 


157.—ANALYTICAL APPROXIMATIONS. [See also Note 153.] The first 
twelve approximations concern the functions e~*Jo(x) and e~*J,(x) in which 
I(x) and J,(x) are the usual Bessel functions of imaginary argument. 


To better than .0006 over (0, ~), 


( 1 + .297x + .341x2 ) 
e*Io(*) 4 2.137x? + 2.0968) 


To better than .00005 over (0, ~), 


1 + .302x + .234x% + ) 
2.2979x + 2.3871x* + 1.2032x° + .7183x4) 


To better than .000,009 over (0, 1), 

e-*Io(x) = (1 + .0302x + .0889x*)/(1 + 1.0299x + .3728x?). 
To better than .0001 over (0, 2), 

e*Ig(x) = (1 + .1693x% + .0844x*)/(1 + 1.1665x + .5247x*). 
To better than .0007 over (0, 4), 

e-*Ig(x) = (1 + .4537x + .0955x?)/(1 + 1.4387x + .8855x*). 
To better than .0005 over (0, ~), 

= «(3.78 + 9.81x + 3.09x? + 

To better than .0003 over (1, ©), 

= «(4.51 + 9.00x + 3.25x? + 

To better than .00005 over (2, ~), 

= «(10.281 + 3.752x + 4.541%? + 6.296x*)-1. 

To better than .000,006 over (0, 1), 

= (.49974x — .01695x*)/(1 + .95935x + .36282x*). 
To better than .00007 over (0, 2), 

(x) = (.4981x + .0066x?)/(1 + .9805x + .4477x*). 

To better than .0003 over (0, 4), 

e-*I,(x) = (.4935x + .0268x*)/(1 + .9667x + .5373x*). 

To better than .00008 over O< x< ~, 

(x) = «(15.4 + 74.8% + 67.2x? + 235.8x° + 43.5x4 


+ 59.4x5 + 39.6x8)-1. 


Definite Integral: To better than .00055 over O< x< ~, 


N(x) = 30x-4 f — dt 
0 


= (1 + .38382x® — .55605x7 + .34791x* — .10369x° + .01245x")—'. 


(48) 
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(52) 


(53) 


(54) 


(55) 
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Mach Number in Terms of Pressure Ratio: To better than .0011 over 
.3< M< 1.0, the inverse of (a) is (b). 
(a) x = Ps/Pe = [1 + (y = 14. 
(b) M = (2.714 — 2.625x)/(1 + 1.650x — 1.955x?). 
Mach Number in Terms of Pressure Ratio: To better than .0014 over 
1< M< 3, the inverse of (c) is (d). 
(c) x = Ps/Pe = (2yM*/(y + 1) — (vy — 1)/(v + 

xX + = 14. 
(d) M = (8.19 + 29.40x — 24.58x*)/(1 + 30x). 


Mach Number in Terms of Pressure Ratio: To better than .0021 
over .3< M< 3, the inverse of (a) over .3< M< 1 conjoined with 
(c) over 1< M< 3 is given by 


M = (8.11 + 23.60x — 39.66x? + 8.98x*)/ 
(1 + 28.70x — 15.99x* — 5.74x*). 


Essentially the same approximation was reported to us by Mr. 
Rapp. 


Pearson Cosine Transformation : To better than .0014 over 0< x< 1, 
r(x) = cos x/(1 + x#) = (-1 — 7.47x + 8.47x*)/ 
(1 + 11.65x + 12.05x*). 
r(x!) = —r(x) can be used to obtain function values over 1 < x< @. 
Pearson Cosine Transformation : To better than .00017 over 0< x< 1, 
—1 — 4.828n + 7.866? — 
1 + 5.560n — 4.985y? + .385y° + .84x 


Natural Addition Logarithm : To better than .00026 overO< x< ~, 
In (1 + e-*) = (1 + .3581x + .1151x? + .0094x° + .0052x*)- In 2. 
Natural Addition Logarithm : To better than .000,045 overO< x< ~, 
In (1 + e-*) = (1 + .36123x + .10204x? + .02411x* 

— .00055x* + .00069x5)-? In 2. 
Natural Addition Logarithm : To better than .000,008 overO< x< @, 


use constants .360571, .105546, .018760, .002654, —.000100, .000066 
in next form of sequence. 


r(x) = 


HASTINGs, JR. 
James P. Wong, Jr. 


RAND Corporation 
1700 Main Street 
Santa Monica, California 


158.—CounTs OF TWIN PRIMES LESS THAN 100000. The writer has 
constructed. a table making the counting of twin primes relatively easy. 
The counts differ from those of previous observers to such an extent that it 
seems desirable to publish the results. 


| 
: 


48 NOTES 


The previous tables are identified merely by giving the author’s names. 
The exact references are given in LEHMER’s Guide. 

Table I shows the various counts of twin primes for each of the 10 
myriads. Following GLAISHER 1 and 3 are counted as pair of twin primes. 


TABLE I 
Myriad SUTTON STACKEL 

1 206 206 
2 137 137 
3 125 125 
4 124 124 
5 113 114 
6 106 106 
7 93 94 
8 102 102 
9 108 109 
10 108 108 
Totals 1222 1225 


The results of GLAISHER agree exactly with those of the author. POLETTI 


POLETTI 


203 


GLAISHER & 
SEXTON 


206 


1225 


gives a total of 1217 but the figures he gives add in fact to 1216. 
The total 1225 is in agreement with Harpy & LiTrLEwoop. 


Table II shows the distribution of these twin primes as to their final 


digits. 
TABLE II 
Digital Endings 
Myriad (1, 3) (7, 9) 

1 68 64 

2 49 46 

3 42 45 

4 46 39 

5 36 40 

6 39 32 

24 36 

8 34 33 

9 33 35 

10 31 44 

Totals 402 414 

POLETTI 401 411 


(9, 1) 


The reference here is to PoLETTI 2, Table XVII. 


The author has made a count of double prime pairs like 191, 193, 197, | 
199. The number 2, of these in the r-th myriad is as follows: 


ee 2 3 4 5 
%: 11 7 3 2 1 


6 
2 


7 
3 


8 
3 


9 
3 


10 
2 


There are in all 37 such double pairs. The set 1, 3, 7, 9 is not counted. 
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Primes of the forms 6n — 1 and 6n + 1 have also been counted and are 
compared with Poletti in Table III. 


TABLE III 
6n -1 6n + 1 
Myriad POLETTI SEXTON POLETTI SEXTON 

1 616 616 614 611 
2 519 520 513 513 
3 498 497 486 486 
4 479 479 479 479 
5 463 463 467 467 
6 466 466 458 458 
7 449 449 429 429 
8 447 447 455 455 
9 437 436 442 440 
10 432 433 444 446 
Totals 4806 4806 4787 4784 


The total 4784 agrees with the count given in CUNNINGHAM 28. 


C. R. SEXTON 
1630 Live Oak Way 
Walnut Creek, California 


1D. H. LenMer, Guide to Tables in the Theory of Numbers. Nat. Res. Council, Bull. 
no. 105. Washington, 1941. 


159.—A NEw APPROXIMATION TO THE RECIPROCAL OF z. In 1949 I de- 
cided to extend to about 1120D my previously published! approximation to 
=, To attain a completely reliable result it was decided to evaluate the 
reciprocal of « by two methods: first, the summation of a suitable highly 
convergent series; and second, machine division of unity by an approxima- 
tion to x that Levi B. Smita and I had computed to about 1120D? in exten- 
sion of an 808D value we had previously published® in collaboration with 
D. F. FERGuson. 

Smith computed by desk machine the sum of the following series due to 
RAMANUJAN :4 


1123 22583 11-3 44043 1-31-3-5-7 
882 2 4 8825 2-4 47-82 


Each term, omitting the coefficients 1123, 22583, ---, in arithmetic pro- 
gression, was evaluated to 1130D and then these coefficients were introduced 
in the calculation to complete the summation. Two independent congru- 
ential checks were applied at intervals of 100D at every stage of this 
calculation. 

Concurrently I used the ENIAC approximation’ to as the divisor in 
an independent evaluation of 2—'. Smith’s final result was compared with 
mine in October 1952 and complete agreement through 1123D was observed, 
thereby verifying to a comparable degree of accuracy the ENIAC deter- 
mination of 7. 


1es. 
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NOTES 


Subsequently I continued the desk machine reciprocation of x to 1408D. 
As an exercise on the UNIVAC the complete product of this approximation 
to x and the ENIAC 2040D work sheet value of + was computed on 20 
April 1953 by two independent routines in an average time of about 8 
minutes. The accuracy of my desk machine calculation was confirmed and 
the extension of x~! to 2037D was at once deduced on the UNIVAC. The 


final result of these investigations is reproduced below. 


x = 0.31830 98861 
28974 95334 
12145 68545 
39967 84584 
57924 34795 
11431 88092 


83790 
68811 
35159 
79933 
32338 
37801 
86800 
60939 
99677 
93095 
65431 
12120 
32481 
78126 
55197 


67153 
77935 
16073 
87481 
67247 
73805 
49734 
88973 
36198 
07397 
35510 
33606 
35512 
61360 
14178 
12291 
32136 
81627 
69869 
52625 
21672 
12622 
90196 
41999 
97859 
74755 
13628 
47799 
11554 
99207 


77675 
95268 
78582 
81551 
80483 
34791 
46954 
04375 
39992 


26745 
45307 
36922 
46155 
44725 
22409 
78919 
76314 
48857 
78920 


02872 
01802 
29157 
49279 
80236 
78821 
22179 
95731 
51170 
72951 
56327 
23974 
62176 
35226 
70263 
17849 
15149 
16509 
61743 
54359 
95162 
28413 
99094 
32258 
15523 
44543 
56557 
49855 
89140 
88007 
36731 
96299 
39373 
58707 


19291 
53250 
59348 
15377 
22844 
56881 
93566 
48207 
77168 
36118 
04589 
55733 
98752 
62813 
00326 
82870 
44519 
66574 
28073 
42251 
88324 
20258 
86324 
66136 
14682 
73251 
68535 
47188 
53958 
97494 
09559 
93116 
89313 
55658 
92036 
05349 
63953 
07900 
58471 
69486 


J. W. WRENCH, JR. 


|_| 
1 
a arcco 
p. 16 
, 3 
p. 24 
4 
Math 
p. 23 
5 
40689 48091 
| 27605 61719 
30575 21463 
4 38506 43478 
64760 53095 | 
87387 71057} 
44619 98928 66461 14981 
23339 72925 33928 79917 
48278 69721 34235 62235 
03753 43210 60194 86675 60498 
89932 70610 06440 64955 94332 34962 
" 39196 33168 07199 62678 99766 08870 Olym 
55951 01400 87776 99142 02443 29536 
27555 29475 92915 95696 24854 99215 Holla 
50049 00059 11380 55935 05042 35492 
04184 96232 12481 24062 92968 69183 42315 
08151 12401 74305 04434 31828 49165 54925 
- 70799 75031 06587 96354 48187 59414 38081 V. 7, 
lg 39995 18153 15415 40787 17965 46851 37902 
. 33250 91411 88665 37300 05224 42306 99008 
77335 89007 52511 63423 39051 56449 66686 
29021 22470 73757 72733 84334 94939 50115 | 
s 66721 06239 21718 79113 43741 93020 76310 
fs 35161 67888 59599 01050 87751 89176 92101 
ee 57058 30302 82080 77012 77632 93986 07799 
ee 91573 83781 19618 44123 75086 78602 05224 
ae 77560 77507 77622 13530 86816 05386 99112 
14158 07721 20705 24902 51991 25940 19116 | 
86023 29659 28237 24811 50889 43579 48189 , 
> 80654 58954 04332 13063 63070 68137 35383 
: 17752 63819 33013 92880 95539 41375 35562 59090 
. 07067 91516 60376 36773 75875 53224 06119 04381 
- 67197 50207 02542 58086 46316 09974 75551 26924 | 
ae 42068 40888 17109 95700 75854 77388 32387 57471 | 
" 87568 69406 46047 42916 75847 11423 72726 83858 63645 
< 83928 33001 75661 58662 70699 55819 94917 29858 01219 
; 78737 81891 76610 06740 61076 10946 24643 16188 52064 | 
: 56626 28379 61949 96448 76670 34871 39796 95002 13677 
60079 57344 71992 16048 00547 80217 49909 70957 36522 
| 27989 78065 37994 85416 69922 29841 65780 75535 07100 
91369 12167 34295 86169 13446 65407 09707 85 
4505 Strathmore Ave. 
Garrett Park, Md. . 
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QUERY 
43. INCOMPLETE HANKEL FUNCTION. Has the Incomplete Hankel function 


f(x,y) = deexp /(@ 
been tabulated? If so, where can such a table be obtained? 


A. P. BurGer 
Olympiaplein 32, 
Amsterdam-Z. 
Holland 


CORRIGENDUM 
V. 7, p. 275, 1. —2, for MaG read Mac. 


D. 
on 
20 ‘ 
8 
nd 
he 
91 
119 
163 
95 
57 
81 
17 
35 
198 
62 
870 
~ 
215 
492 
315 
081 
902 
008 
686 
15 | 
310 
101 
799 
224 
112 
116 | 
189 ' 
383 
090 
381 - 
924 | 
471 
645 
219 
064 — 
677 
522 
1100 
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